在本文中，将探讨高斯过程回归（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"初始:\n{kernel}\n最优:\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"初始:\n{kernel}\n最优:\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)
正如所希望的，随机重启允许优化在不良初始值的情况下找到最佳的超参数集。