I am using sklearn gaussian processes to fit my data, and am trying to sample the posterior in order to do some uncertainty estimation. However, it seems that using the function .sample_y is sampling the prior, because although the gaussian process fit to my data is good, the samples are all over the place..
I define my kernels and then initiate my GP regressor, and fit to the data:
gp_regressor = GaussianProcessRegressor(kernel=kernel, alpha=data_yerr_scaled_fit**2, random_state=0,n_restarts_optimizer=n_optimizer).fit(stacked_data_x_for_fit, data_y_scaled_fit)
I create my best fitting model by using the predict function:
gp, gp_err = gp_regressor.predict(stacked_data, return_std=True)
and then I try to sample the regressor (which I thought had been trained at this point) but this returns samples of the prior:
samples=gp_regressor.sample_y(stacked_data_x,n_samples=number_samples, random_state=0)
How can I make sure I am sampling the posterior?