在 Python 中再现 R 的高斯过程最大似然回归
Reproducing R's gaussian process maximum likelihood regression in Python
我在 R 中实现了一个函数来估计基本 sin 函数的高斯过程参数。不幸的是，该项目必须在 Python 中完成，我一直在尝试使用 SKlearn 在 python 中重现 R 库的 hetGP 的行为，但我很难将前者映射到后者。
我对高斯过程的理解仍然有限，而且我是 sklearn 的初学者，所以我非常感谢对此的一些帮助。
我的 R 代码：
library(hetGP)
set.seed(123)
nvar <- 2
n <- 400
r <- 1
f <- function(x) sin(sum(x))
true_C <- matrix(1/8 * (3 + 2 * cos(2) - cos(4)), nrow = 2, ncol = 2)
design <- matrix(runif(nvar*n), ncol = nvar)
response <- apply(design, 1, f)
model <- mleHomGP(design, response, lower = rep(1e-4, nvar), upper = rep(1,nvar))
稍后在代码中，我使用 model$Ki 和 model$theta
model$theta: 0.9396363 0.9669170
dim(model$ki): 400 400
到目前为止我的 Python 代码：
import numpy as np
from sklearn.gaussian_process import GaussianProcessRegressor
from sklearn.gaussian_process.kernels import RBF
n = 400
n_var = 2
real_c = np.full((2, 2), 1 / 8 * (3 + 2 * np.cos(2) - np.cos(4)))
design = np.random.uniform(size=n * n_var).reshape(-1, 2)
test = np.random.uniform(size=n * n_var).reshape(-1, 2)
response = np.apply_along_axis(lambda x: np.sin(np.sum(x)), 1, design)
kernel = RBF(length_scale=(1, 1))
gpr = GaussianProcessRegressor(kernel=kernel, n_restarts_optimizer=10,
optimizer="fmin_l_bfgs_b").fit(design, response)
gpr.predict(test, return_std=True)
theta = gpr.kernel_.get_params()["length_scale"]
#theta = gpr.kernel_.theta
k_inv = gpr._K_inv
theta = [1.78106558 1.80083585]
经过一个多星期，我终于找到了答案。
（通过查看 Scikit-learn 和 hetGP 实现）
他们的实现有很多不同点。
- 首先，必须在 sklearn 中显式实例化噪声级别和西格玛，否则它们将不会被优化。
- 此外，找到k_inv的最佳方法是使用GaussianProcessRegressor.L_，K的Cholesky分解的下部
- 最后
hetGP 没有按 sigma 缩放他们的 K，所以我们必须手动完成。
我是怎么做的：
import scipy
import numpy as np
from sklearn.gaussian_process import GaussianProcessRegressor
from sklearn.gaussian_process.kernels import ConstantKernel, RBF, WhiteKernel
n = 400
n_var = 2
design = np.random.uniform(size=n * n_var).reshape((-1, n_var))
response = np.array([np.sin(np.sum(row)) for row in design])
real_c = np.full((n_var, n_var), 1 / 8 * (3 + 2 * np.cos(2) - np.cos(4)))
# Kernel has to have a constant kernel for sigma, a n_var dimension length_scale and a white kernel to optimize the noise
kernel = ConstantKernel(1e-8) * RBF(length_scale=np.array([1.0, 1.0])) + WhiteKernel(noise_level=1)
gpr = GaussianProcessRegressor(kernel=kernel, alpha=1e-10)
gpr.fit(design, response)
L_inv = scipy.linalg.solve_triangular(gpr.L_.T, np.eye(gpr.L_.shape[0]))
k_inv = L_inv.dot(L_inv.T)
sigma_f = gpr.kernel_.k1.get_params()['k1__constant_value']
theta = gpr.kernel_.k1.get_params()['k2__length_scale']
Ki = k_inv * sigma_f
我在 R 中实现了一个函数来估计基本 sin 函数的高斯过程参数。不幸的是，该项目必须在 Python 中完成，我一直在尝试使用 SKlearn 在 python 中重现 R 库的 hetGP 的行为，但我很难将前者映射到后者。
我对高斯过程的理解仍然有限，而且我是 sklearn 的初学者，所以我非常感谢对此的一些帮助。
我的 R 代码：
library(hetGP)
set.seed(123)
nvar <- 2
n <- 400
r <- 1
f <- function(x) sin(sum(x))
true_C <- matrix(1/8 * (3 + 2 * cos(2) - cos(4)), nrow = 2, ncol = 2)
design <- matrix(runif(nvar*n), ncol = nvar)
response <- apply(design, 1, f)
model <- mleHomGP(design, response, lower = rep(1e-4, nvar), upper = rep(1,nvar))
稍后在代码中，我使用 model$Ki 和 model$theta
model$theta: 0.9396363 0.9669170
dim(model$ki): 400 400
到目前为止我的 Python 代码：
import numpy as np
from sklearn.gaussian_process import GaussianProcessRegressor
from sklearn.gaussian_process.kernels import RBF
n = 400
n_var = 2
real_c = np.full((2, 2), 1 / 8 * (3 + 2 * np.cos(2) - np.cos(4)))
design = np.random.uniform(size=n * n_var).reshape(-1, 2)
test = np.random.uniform(size=n * n_var).reshape(-1, 2)
response = np.apply_along_axis(lambda x: np.sin(np.sum(x)), 1, design)
kernel = RBF(length_scale=(1, 1))
gpr = GaussianProcessRegressor(kernel=kernel, n_restarts_optimizer=10,
optimizer="fmin_l_bfgs_b").fit(design, response)
gpr.predict(test, return_std=True)
theta = gpr.kernel_.get_params()["length_scale"]
#theta = gpr.kernel_.theta
k_inv = gpr._K_inv
theta = [1.78106558 1.80083585]
经过一个多星期，我终于找到了答案。
（通过查看 Scikit-learn 和 hetGP 实现）
他们的实现有很多不同点。
- 首先，必须在 sklearn 中显式实例化噪声级别和西格玛，否则它们将不会被优化。
- 此外，找到k_inv的最佳方法是使用GaussianProcessRegressor.L_，K的Cholesky分解的下部
- 最后
hetGP没有按sigma缩放他们的 K，所以我们必须手动完成。
我是怎么做的：
import scipy
import numpy as np
from sklearn.gaussian_process import GaussianProcessRegressor
from sklearn.gaussian_process.kernels import ConstantKernel, RBF, WhiteKernel
n = 400
n_var = 2
design = np.random.uniform(size=n * n_var).reshape((-1, n_var))
response = np.array([np.sin(np.sum(row)) for row in design])
real_c = np.full((n_var, n_var), 1 / 8 * (3 + 2 * np.cos(2) - np.cos(4)))
# Kernel has to have a constant kernel for sigma, a n_var dimension length_scale and a white kernel to optimize the noise
kernel = ConstantKernel(1e-8) * RBF(length_scale=np.array([1.0, 1.0])) + WhiteKernel(noise_level=1)
gpr = GaussianProcessRegressor(kernel=kernel, alpha=1e-10)
gpr.fit(design, response)
L_inv = scipy.linalg.solve_triangular(gpr.L_.T, np.eye(gpr.L_.shape[0]))
k_inv = L_inv.dot(L_inv.T)
sigma_f = gpr.kernel_.k1.get_params()['k1__constant_value']
theta = gpr.kernel_.k1.get_params()['k2__length_scale']
Ki = k_inv * sigma_f