SKlearn Gaussian Process with constant, manually set correlation

Question

I want to use the Gaussian Process approximation for a simple 1D test function to illustrate a few things. I want to iterate over a few different values for the correlation matrix (since this is 1D it is just a single value) and show what effect different values have on the approximation. My understanding is, that "theta" is the parameter for this. Therefore I want to set the theta value manually and don't want any optimization/changes to it. I thought the constant kernel and the clone_with_theta function might get me what I want but I didn't get it to work. Here is what I have so far:

import numpy as np
from sklearn.gaussian_process import GaussianProcessRegressor
from sklearn.gaussian_process.kernels import RBF, ConstantKernel as ConstantKernel


def f(x):
"""The function to predict."""
return x/2 + ((1/10 + x)  * np.sin(5*x - 1))/(1 + x**2 * (np.sin(x - (1/2))**2))

# ----------------------------------------------------------------------
#  Data Points
X = np.atleast_2d(np.delete(np.linspace(-1,1, 7),4)).T
y = f(X).ravel()

# Instantiate a Gaussian Process model
kernel = ConstantKernel(constant_value=1, constant_value_bounds='fixed')

theta = np.array([0.5,0.5])
kernel = kernel.clone_with_theta(theta)

gp = GaussianProcessRegressor(kernel=kernel, optimizer=None)

# Fit to data using Maximum Likelihood Estimation of the parameters
gp.fit(X, y)

# Make the prediction on the meshed x-axis (ask for MSE as well)
y_pred, sigma = gp.predict(x, return_std=True)

# Plot 
# ...

Answer 1

I programmed a simple implementation myself now, which allows to set correlation (here 'b') manually:

import numpy as np
from numpy.linalg import inv
    
def f(x):
    """The function to predict."""
    return x/2 + ((1/10 + x)  * np.sin(5*x - 1))/(1 + x**2 * (np.sin(x - (1/2))**2))

def kriging_approx(x,xt,yt,b,mu,R_inv):
    
    N = yt.size
    one = np.matrix(np.ones((yt.size))).T
    
    r = np.zeros((N))
    for i in range(0,N):
        r[i]= np.exp(-b * (xt[i]-x)**2)
    

    y = mu + np.matmul(np.matmul(r.T,R_inv),yt - mu*one)
    y = y[0,0]
    return y

def calc_R (x,b):
    N = x.size
    # setup R
    R = np.zeros((N,N))
    for i in range(0,N):
        for j in range(0,N):
            R[i][j] = np.exp(-b * (x[i]-x[j])**2)
            
    R_inv = inv(R)
    return R, R_inv

def calc_mu_sig (yt, R_inv):
    N = yt.size
    one = np.matrix(np.ones((N))).T
    mu = np.matmul(np.matmul(one.T,R_inv),yt) / np.matmul(np.matmul(one.T,R_inv),one)
    mu = mu[0,0]
    
    sig2 = (np.matmul(np.matmul((yt - mu*one).T,R_inv),yt - mu*one))/(N)
    sig2 = sig2[0,0]
    
    return mu, sig2

# ----------------------------------------------------------------------

#  Data Points
xt = np.linspace(-1,1, 7)
yt = np.matrix((f(xt))).T

# Calc R   
R, R_inv = calc_R(xt, b)

# Calc mu and sigma   
mu_dach, sig_dach2 = calc_mu_sig(yt, R_inv)

# Point to get approximation for
x = 1

y_approx = kriging_approx(x, xt, yt, b, mu_dach, R_inv)