如何为 scikit-learn 的高斯过程回归指定先验？

Posted 2023-03-12

【中文标题】如何为 scikit-learn 的高斯过程回归指定先验？

【英文标题】：How to specify the prior for scikit-learn's Gaussian process regression?

【发布时间】：2019-07-23 21:44:16

【问题描述】：

如here 所述，scikit-learn 的高斯过程回归 (GPR) 允许“无需预先拟合的预测（基于 GP 先验）”。但我对我的先验应该是什么有一个想法（即它不应该简单地具有零平均值，但也许我的输出 y 与我的输入 X 成线性比例，即 y = X）。我如何将这些信息编码到 GPR 中？

下面是一个工作示例，但它假设我之前的均值为零。我read 说“GaussianProcessRegressor 不允许指定均值函数，总是假设它是零函数，突出了均值函数在计算后验中的作用减弱。”我相信这是custom kernels（例如异方差）背后的动机，在不同的X 上具有可变比例，尽管我仍在努力更好地了解它们提供的功能。有没有办法绕过零均值先验，以便可以在 scikit-learn 中指定任意先验？

import numpy as np
from matplotlib import pyplot as plt
from sklearn.gaussian_process import GaussianProcessRegressor
from sklearn.gaussian_process.kernels import RBF, ConstantKernel as C

def f(x):
    """The function to predict."""
    return 1.5*(1. - np.tanh(100.*(x-0.96))) + 1.5*x*(x-0.95) + 0.4 + 1.5*(1.-x)* np.random.random(x.shape)

# Instantiate a Gaussian Process model
kernel = C(10.0, (1e-5, 1e5)) * RBF(10.0, (1e-5, 1e5))

X = np.array([0.803,0.827,0.861,0.875,0.892,0.905,
                0.91,0.92,0.925,0.935,0.941,0.947,0.96,
                0.974,0.985,0.995,1.0])
X = np.atleast_2d(X).T

# Observations and noise
y = f(X).ravel() 
noise = np.linspace(0.4,0.3,len(X))
y += noise

# Instantiate a Gaussian Process model
gp = GaussianProcessRegressor(kernel=kernel, alpha=noise ** 2,
                              n_restarts_optimizer=10)
# 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)
x = np.atleast_2d(np.linspace(0.8, 1.02, 1000)).T
y_pred, sigma = gp.predict(x, return_std=True)

plt.figure() 
plt.errorbar(X.ravel(), y, noise, fmt='k.', markersize=10, label=u'Observations')
plt.plot(x, y_pred, 'k-', label=u'Prediction')
plt.fill(np.concatenate([x, x[::-1]]),
         np.concatenate([y_pred - 1.9600 * sigma,
                        (y_pred + 1.9600 * sigma)[::-1]]),
         alpha=.1, fc='k', ec='None', label='95% confidence interval')
plt.xlabel('x')
plt.ylabel('y')
plt.xlim(0.8, 1.02)
plt.ylim(0, 5)
plt.legend(loc='lower left')
plt.show()

【问题讨论】：

如果它可能有任何用处，我从您调用 gp.fit(X, y) 中找到了一个很好的数据近似方程，它是一个标准的四参数逻辑方程，“y = d + (a - d) / (1.0 + pow(X / c, b))" 参数 a = 3.7534001422009142E+00, b = 1.9263507395092211E+02, c = 9.6018856179845813E-01, d = 7.8032402821316 -01 产生 R 平方 = 0.9995 和 RMSE = 0.03

@JamesPhillips 感谢您的意见。我刚刚根据上面定义的f(x) 创建了带有添加噪声的任意数据（这应该类似于您指出的修改后的 tanh），这里使用的数据对我来说并不重要。我只是对能够使用我可以控制的全科医生感兴趣。

【参考方案1】：

这是一个关于如何将先验均值函数用于 sklearn GPR 模型的示例。

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

A=np.linspace(5,25,num=100)
# prior mean function
prior_beta=12-0.3*A
# true function
true_beta=20-0.7*A

rng = np.random.seed(44)
# Training data
size=15
ind=np.random.randint(0,100,size=size)
# generate the posterior variance (noisy samples)
var_=np.random.uniform(0.1,10.0,size=size)
A_=A[ind][:, np.newaxis]
beta_=true_beta[ind]-prior_beta[ind]
beta_1=true_beta[ind]

plt.figure()

kernel = ConstantKernel(4) * RBF(length_scale=2, length_scale_bounds=(1e-3, 1e2))
gp = GaussianProcessRegressor(kernel=kernel,
                              alpha=var_,optimizer=None).fit(A_, beta_)
X_ = np.linspace(5, 25, 100)
y_mean, y_cov = gp.predict(X_[:, np.newaxis], return_cov=True)
# Now you add the prior mean function back
y_mean=y_mean+12-0.3*X_
plt.plot(X_, y_mean, 'k', lw=3, zorder=9, label='predicted')
plt.fill_between(X_, y_mean - 3*np.sqrt(np.diag(y_cov)),
                 y_mean + 3*np.sqrt(np.diag(y_cov)),
                 alpha=0.5, color='k', label='+-3sigma')
plt.plot(A,true_beta, 'r', lw=3, zorder=9,label='truth')
plt.plot(A,prior_beta, 'blue', lw=3, zorder=9,label='prior')
plt.errorbar(A_[:,0], beta_1, yerr=3*np.sqrt(var_), fmt='x',ecolor='g',marker='s', 
mfc='g', ms=10,capsize=6,label='training set')

plt.title("Initial: %s\n"% (kernel))
plt.legend()
plt.show()