Calculate the probability of a function given by a Gaussian process model

I am using Gaussian process regression using scikit-learn. (my actually simple one-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 was fit, I would like to calculate the likelihood that another dataset ( y1

) will be under this installed model. eg P(y1|GP_fit_on_y)

. It would be very naive.

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 GP covariance structure. Any idea on how to do this?