高斯过程回归（GPR）是一种强大的机器学习技术，它不仅能够预测连续值，还能提供预测的不确定性。本文将探讨GPR在估计数据中噪声水平的能力，并展示核函数超参数初始化的重要性。
将在一个单特征的环境中工作，创建一个目标生成函数，该函数将生成待预测的目标，并提供一个添加噪声的选项。
import numpy as np
def target_generator(X, add_noise=False):
target = 0.5 + np.sin(3 * X)
if add_noise:
rng = np.random.RandomState(1)
target += rng.normal(0, 0.3, size=target.shape)
return target.squeeze()
首先，观察不添加任何噪声的目标生成器，以观察希望预测的信号。
X = np.linspace(0, 5, num=80).reshape(-1, 1)
y = target_generator(X, add_noise=False)
import matplotlib.pyplot as plt
plt.plot(X, y, label="预期信号")
plt.legend()
plt.xlabel("X")
plt.ylabel("y")
目标函数通过正弦函数转换输入X。现在，将生成一些带噪声的训练样本。为了说明噪声水平，将真实信号与带噪声的训练样本一起绘制。
rng = np.random.RandomState(0)
X_train = rng.uniform(0, 5, size=20).reshape(-1, 1)
y_train = target_generator(X_train, add_noise=True)
plt.plot(X, y, label="预期信号")
plt.scatter(x=X_train[:, 0], y=y_train, color="black", alpha=0.4, label="观测值")
plt.legend()
plt.xlabel("X")
plt.ylabel("y")
接下来，将创建一个使用加性核的高斯过程回归模型，该模型结合了径向基函数（RBF）核和白噪声核（WhiteKernel）。白噪声核能够估计数据中存在的噪声量，而RBF核则用于拟合数据和目标之间的非线性关系。然而，将展示超参数空间包含多个局部最小值，这突出了初始超参数值的重要性。
from sklearn.gaussian_process import GaussianProcessRegressor
from sklearn.gaussian_process.kernels import RBF, WhiteKernel
kernel = 1.0 * RBF(length_scale=1e1, length_scale_bounds=(1e-2, 1e3)) + WhiteKernel(noise_level=1, noise_level_bounds=(1e-10, 1e1))
gpr = GaussianProcessRegressor(kernel=kernel, alpha=0.0)
gpr.fit(X_train, y_train)
y_mean, y_std = gpr.predict(X, return_std=True)
观察到，尽管找到了最优核，但噪声水平仍然很高，长度尺度甚至更大。长度尺度达到了为该参数允许的最大范围，因此收到了一个警告。更重要的是，注意到模型没有提供有用的预测：平均预测似乎是恒定的，它没有遵循预期的无噪声信号。
plt.plot(X, y, label="预期信号")
plt.scatter(x=X_train[:, 0], y=y_train, color="black", alpha=0.4, label="观测值")
plt.errorbar(X, y_mean, y_std, label="后验均值 ± 标准差")
plt.legend()
plt.xlabel("X")
plt.ylabel("y")
plt.title(f"初始: {kernel}\n最优: {gpr.kernel_}\n对数边际似然: {gpr.log_marginal_likelihood(gpr.kernel_.theta)}", fontsize=8)
现在，将RBF的初始长度尺度值设得更大，白噪声核的初始噪声水平更低，同时保持参数范围不变。
kernel = 1.0 * RBF(length_scale=1e-1, length_scale_bounds=(1e-2, 1e3)) + WhiteKernel(noise_level=1e-2, noise_level_bounds=(1e-10, 1e1))
gpr = GaussianProcessRegressor(kernel=kernel, alpha=0.0)
gpr.fit(X_train, y_train)
y_mean, y_std = gpr.predict(X, return_std=True)
首先，观察到模型的预测比之前的模型更精确：这个新模型能够估计无噪声的函数关系。查看核超参数，发现找到的最佳组合比第一个模型具有更小的噪声水平和更短的长度尺度。
plt.plot(X, y, label="预期信号")
plt.scatter(x=X_train[:, 0], y=y_train, color="black", alpha=0.4, label="观测值")
plt.errorbar(X, y_mean, y_std, label="后验均值 ± 标准差")
plt.legend()
plt.xlabel("X")
plt.ylabel("y")
plt.title(f"初始: {kernel}\n最优: {gpr.kernel_}\n对数边际似然: {gpr.log_marginal_likelihood(gpr.kernel_.theta)}", fontsize=8)
可以为不同的超参数检查高斯过程回归的对数边际似然（LML），以了解局部最小值。
from matplotlib.colors import LogNorm
length_scale = np.logspace(-2, 4, num=80)
noise_level = np.logspace(-2, 1, num=80)
length_scale_grid, noise_level_grid = np.meshgrid(length_scale, noise_level)
log_marginal_likelihood = [gpr.log_marginal_likelihood(theta=np.log([0.36, scale, noise])) for scale, noise in zip(length_scale_grid.ravel(), noise_level_grid.ravel())]
log_marginal_likelihood = np.reshape(log_marginal_likelihood, noise_level_grid.shape)
vmin, vmax = (-log_marginal_likelihood).min(), 50
level = np.around(np.logspace(np.log10(vmin), np.log10(vmax), num=20), decimals=1)
plt.contour(length_scale_grid, noise_level_grid, -log_marginal_likelihood, levels=level, norm=LogNorm(vmin=vmin, vmax=vmax))
plt.colorbar()
plt.xscale("log")
plt.yscale("log")
plt.xlabel("长度尺度")
plt.ylabel("噪声水平")
plt.title("对数边际似然")
plt.show()
可以看到，有两个局部最小值对应于之前找到的超参数组合。根据超参数的初始值，基于梯度的优化可能会或可能不会收敛到最佳模型。因此，重复多次优化对于不同的初始化是很重要的。这可以通过设置高斯过程回归类的n_restarts_optimizer参数来完成。
kernel = 1.0 * RBF(length_scale=1e1, length_scale_bounds=(1e-2, 1e3)) + WhiteKernel(noise_level=1, noise_level_bounds=(1e-10, 1e1))
gpr = GaussianProcessRegressor(kernel=kernel, alpha=0.0, n_restarts_optimizer=10, random_state=0)
gpr.fit(X_train, y_train)
y_mean, y_std = gpr.predict(X, return_std=True)
正如所希望的，随机重启允许优化在不良初始值的情况下找到最佳超参数集合。