具有恒定的手动设置相关性的 SKlearn 高斯过程答案

【问题标题】：SKlearn Gaussian Process with constant, manually set correlation具有恒定的手动设置相关性的 SKlearn 高斯过程
【发布时间】：2020-11-21 17:17:35
【问题描述】：

我想对简单的一维测试函数使用高斯过程近似来说明一些事情。我想迭代相关矩阵的几个不同值（因为这是一维的，它只是一个值）并显示不同值对近似值的影响。我的理解是，“theta”是这个参数。因此，我想手动设置 theta 值，并且不希望对其进行任何优化/更改。我认为常量内核和 clone_with_theta 函数可能会得到我想要的东西，但我没有让它工作。这是我目前所拥有的：

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 
# ...

【问题讨论】：

标签： python scikit-learn gaussian-process

【解决方案1】：

我现在自己编写了一个简单的实现，它允许手动设置相关性（此处为“b”）：

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)

【讨论】：