SKlearn高斯过程具有常数,手动设置相关性



我想对一个简单的1D测试函数使用高斯过程近似来说明一些事情。我想迭代相关矩阵的几个不同值(因为这是1D,所以它只是一个值(,并显示不同值对近似值的影响。我的理解是;θ";是此的参数。因此,我想手动设置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 
# ...

我现在自己编程了一个简单的实现,它允许手动设置相关性(此处为"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)

最新更新