Implementing Cross-Validation for Gaussian Process Regression Although Gaussian Process Module in sklearn package offers an "automatic" optimization based on the posterior likelihood function, I'd like to use cross-validation to pick the best hyperparameters for GP regression model. Now, I met one confusion when using GridSearchCV. Here are two versions of my cross-validation for GP regression (I wrote an auxiliary function "cross_val_kernel" to help explain my code but this is not the key point):
# LOO strategy Cross-Validation, not an efficient algorithm, just for illustration here
def cross_val_kernel(kernels, X, Y):
    performance = {}
    
    for kernel in kernels: 
        likelihood = 0
        for i in range(Y.size):
            gp = GaussianProcessRegressor(kernel = kernel)
            X_train = np.delete(X, i, axis=0)
            Y_train = np.delete(Y, i, axis=0)
            gp.fit(X_train, Y_train)
            y_mean, y_std = gp.predict(X[[i], :], return_std = True)
            likelihood += -np.log(y_std[0]) - (Y[i] - y_mean[0])**2 / (2 * y_std[0]**2)
            
        performance[likelihood] = kernel
        
    return performance

Version 1:
parRange = [np.arange(1, 100, 1), np.arange(0.1, 2, 0.1), np.arange(0.1, 1, 0.1)]
kernels = [ConstantKernel(a, constant_value_bounds='fixed') * RBF(b, length_scale_bounds='fixed') 
           + WhiteKernel(c, noise_level_bounds='fixed') for a, b, c in list(itertools.product(*parRange))]

cross_val_kernel(kernels, X, Y)

Version 2:
parRange = [np.arange(1, 100, 1), np.arange(0.1, 2, 0.1), np.arange(0.1, 1, 0.1)]
kernels = [ConstantKernel(a) * RBF(b) 
           + WhiteKernel(c) for a, b, c in list(itertools.product(*parRange))]

cross_val_kernel(kernels, X, Y)

The only difference between two versions are hyperparameters are fixed or not. I'm wondering which version should I choose for CV, because I found no matter what the parameter I set at first, if the bound is not fixed, the predicted value on the same data set seems to be the same. Thanks!
 A: If the hyperparameters are not fixed, sklearn GP's optimizer will change them to optimal values, so that they best fit the dataset, that's why the prediction obtained is expected to be close to the original values.
For hyper-parameter tuning with grid-search CV, we are explicitly providing the hyper-parameter values and thus will not want sklearn GP optimizer to change the hyper-parameter values, that's why they must be fixed.
Let's try both the versions with Friedman's dataset to make it clearer.
from sklearn.datasets import make_friedman2
X, Y = make_friedman2(n_samples=500, noise=0, random_state=0)

For example, with version 1, as can be seen from the below code, the hyperparameters are not changed by the optimizer and that's what we intend to do if we want explicit hyperpamater tuning.
kernels = [ConstantKernel(a, constant_value_bounds='fixed') * RBF(b, length_scale_bounds='fixed') 
           + WhiteKernel(c, noise_level_bounds='fixed') for a, b, c in list(itertools.product(*parRange))]
gp = GaussianProcessRegressor(kernel = kernels[0])
gp.fit(X, Y)
print(gp.kernel_) 
# 1**2 * RBF(length_scale=0.1) + WhiteKernel(noise_level=0.1)

As opposed to this, with version 2, the GP regressor optimizes the hyperparameter values, that's not something we shall like it to do, since we are computing the marginal-log-likelihoods on our own and want to do model selection based on that. Notice the change in hyperparameter (RBF variance and noise variance) values after fitting with the following code:
kernels = [ConstantKernel(a) * RBF(b) 
           + WhiteKernel(c) for a, b, c in list(itertools.product(*parRange))]
gp = GaussianProcessRegressor(kernel = kernels[0])
gp.fit(X, Y)
print(gp.kernel_)
# 316**2 * RBF(length_scale=0.1) + WhiteKernel(noise_level=1e+05)

