I am fitting a Gaussian process regression using scikit-learn. (mine is actually a simple 1 dimensional case)
from sklearn.gaussian_process import GaussianProcessRegressor
from sklearn.gaussian_process.kernels import RBF, WhiteKernel
import numpy as np
# Some example data
x = np.linspace(0,10,100)
y = np.sin(x) + np.random.normal(0, 0.1, x.shape)
# Fitting the model
some_kernel = RBF() + WhiteKernel()
gpr = GaussianProcessRegressor(kernel=some_kernel)
gpr.fit(x[:, None], y)
Once this model has been fit I would like to calculate the probability that another data set (y1) has under this fitted model. e.g P(y1|GP_fit_on_y).
A very naive way would be.
from scipy.stats import norm as normal
y_pred, std = gpr.predict(x[:, None], return_std=True)
sum_log_lik = np.sum( normal.logpdf(loc=y_pred, scale=std, x=y1) )
But this would not take into account the covariance structure of the GP. Any idea on how to do it?