Python 高斯过程回归:让模型在“不知道”时主动变得不确定
- 2026-08-18 23:11:47
普通回归模型通常只给出一个预测值，但真实分析往往还需要回答另一个问题：这个预测究竟有多可靠？今天用 Python 的 GaussianProcessRegressor 完成一次可复现的高斯过程回归，并故意在训练数据中留出一段空白区，观察模型能否在缺少数据支持的位置扩大不确定性。
今天要解决的问题
我们模拟一个平滑但非线性的连续结局。80 个训练样本只分布在 0–4 和 7–10，中间 4–7 完全没有训练观测；240 个独立测试样本则覆盖整个 0–10 区间。
这不是单纯比较谁的 RMSE 更低，而是同时检查三个问题：
后验均值能否恢复非线性关系； 预测标准差是否会在训练空白区增大； 近似 95% 预测区间能否覆盖相应比例的新观测。
这种思路适合样本量不太大、输入维度较低，并且“预测不确定性”本身具有决策价值的回归任务。
方法原理
高斯过程回归不是预先写死一条直线或某个固定次数的多项式，而是对“可能的函数”设置概率分布。核函数负责描述两个输入位置有多相似：位置越近，函数值通常越相关。
本例使用常数核、RBF 核与白噪声核的组合。RBF 核可写成：
k(x, x′) = σf² exp[−(x − x′)² / (2ℓ²)]
其中，σf 控制函数变化的整体幅度，长度尺度 ℓ 控制相关性随距离衰减的速度，白噪声核则吸收观测噪声。模型拟合后，对每个新位置同时返回后验均值和标准差；在高斯假设下，可用“均值 ± 1.96 × 标准差”构造近似 95% 预测区间。
关键点是：训练数据密集处的函数形状受到更多约束，而远离观测的位置存在更多仍与数据相容的函数，因此后验不确定性通常会增大。
数据与实现思路
完整流程分为三步：先制造一个明确可审计的训练空白区，再通过对数边际似然优化核超参数，最后用独立测试集分别评估有数据支持区和训练空白区。
为了保证结果可复现，代码同时固定 Python、NumPy 和模型内部随机状态。图中的文字全部使用英文，以减少跨平台字体问题；文章解释仍使用中文。
Python 完整实战
构造训练空白区与测试数据
import osfrom pathlib import Path# 把Matplotlib缓存放到可写临时目录，保证代码在受限环境中无警告运行os.environ[”MPLCONFIGDIR”] = ”/tmp/wechat_python_gpr_mpl”os.environ[”MPLBACKEND”] = ”Agg”Path(os.environ[”MPLCONFIGDIR”]).mkdir(parents=True, exist_ok=True)import randomimport numpy as npimport pandas as pdimport matplotlib.pyplot as pltfrom sklearn.gaussian_process import GaussianProcessRegressorfrom sklearn.gaussian_process.kernels import ConstantKernel, RBF, WhiteKernelfrom sklearn.metrics import mean_squared_error, mean_absolute_error, r2_score# 同时固定Python、NumPy和模型内部随机状态SEED = 20260802random.seed(SEED)rng = np.random.default_rng(SEED)# 定义只用于模拟的真实函数，训练点故意避开4到7之间的空白区def true_function(x):return 1.2 * np.sin(1.15 * x) + 0.18 * xx_train = np.concatenate([rng.uniform(0.0, 4.0, 42),rng.uniform(7.0, 10.0, 38),])y_train = true_function(x_train) + rng.normal(0.0, 0.25, x_train.size)# 测试点覆盖完整区间，既评估有训练支持的位置，也评估中间空白区x_test = rng.uniform(0.0, 10.0, 240)y_test = true_function(x_test) + rng.normal(0.0, 0.25, x_test.size)train_gap_count = int(np.sum((x_train > 4.0) & (x_train < 7.0)))test_gap_count = int(np.sum((x_test > 4.0) & (x_test < 7.0)))# 输出样本规模和训练空白区审计，避免把插值与外推混为一谈data_audit = pd.DataFrame({”metric”: [”train_rows”, ”test_rows”, ”train_gap_rows”, ”test_gap_rows”, ”noise_sd”],”value”: [x_train.size, x_test.size, train_gap_count, test_gap_count, 0.25],})print(data_audit.to_string(index=False))
输出：
metric valuetrain_rows 80.00test_rows 240.00train_gap_rows 0.00test_gap_rows 77.00noise_sd 0.25
审计结果确认：训练空白区内没有任何训练样本，但测试集中有 77 个样本落在这里。因此，后续可以直接比较“有训练数据支持”和“没有训练数据支持”两种情形。
拟合高斯过程并绘制后验曲线
# 常数核控制整体幅度，RBF核控制相似性衰减，白噪声核估计观测噪声kernel = (ConstantKernel(1.0, (0.1, 10.0))* RBF(length_scale=1.2, length_scale_bounds=(0.1, 10.0))+ WhiteKernel(noise_level=0.08, noise_level_bounds=(1e-3, 1.0)))gpr = GaussianProcessRegressor(kernel=kernel,normalize_y=True,n_restarts_optimizer=6,random_state=SEED,)# 拟合时通过最大化对数边际似然选择核超参数gpr.fit(x_train.reshape(-1, 1), y_train)learned_amplitude = float(np.sqrt(gpr.kernel_.k1.k1.constant_value))learned_length = float(gpr.kernel_.k1.k2.length_scale)learned_noise_sd = float(np.sqrt(gpr.kernel_.k2.noise_level))print(f”optimized_kernel: {gpr.kernel_}”)print(f”log_marginal_likelihood: {gpr.log_marginal_likelihood_value_:.3f}”)print(f”learned_amplitude: {learned_amplitude:.3f}”)print(f”learned_length_scale: {learned_length:.3f}”)print(f”learned_noise_sd: {learned_noise_sd:.3f}”)# 在规则网格上同时返回后验均值和标准差，构造95%预测带x_grid = np.linspace(0.0, 10.0, 500)grid_mean, grid_std = gpr.predict(x_grid.reshape(-1, 1), return_std=True)grid_lower = grid_mean - 1.96 * grid_stdgrid_upper = grid_mean + 1.96 * grid_std# 第一张图突出训练空白区中的不确定性扩张，而不只展示一条拟合线fig, ax = plt.subplots(figsize=(9.0, 5.5))ax.fill_between(x_grid, grid_lower, grid_upper, color=”#3C78A8”, alpha=0.22,label=”Approx. 95% predictive interval”)ax.plot(x_grid, grid_mean, color=”#28666E”, linewidth=2.2, label=”GPR mean”)ax.plot(x_grid, true_function(x_grid), color=”#B64D46”, linestyle=”--”,linewidth=1.5, label=”Simulation truth”)ax.scatter(x_train, y_train, color=”#263238”, s=28, alpha=0.72,label=”Training observations”, zorder=3)ax.axvspan(4.0, 7.0, color=”#D18B2C”, alpha=0.08, label=”Training gap”)ax.set(title=”Gaussian process regression with location-dependent uncertainty”,xlabel=”Input x”, ylabel=”Outcome y”)ax.legend(loc=”upper left”, frameon=False, ncol=2)ax.grid(alpha=0.20)fig.tight_layout()fig.savefig(”gpr_posterior_curve.png”, dpi=180, bbox_inches=”tight”, facecolor=”white”)plt.close(fig)
输出：
optimized_kernel: 1.23**2 * RBF(length_scale=1.26) + WhiteKernel(noise_level=0.141)log_marginal_likelihood: -50.412learned_amplitude: 1.232learned_length_scale: 1.257learned_noise_sd: 0.375
优化后的长度尺度约为 1.257，说明模型学到的是较平滑、但仍允许局部起伏的函数。橙色背景标出的 4–7 区间没有训练点，蓝色预测带在这里明显变宽；重新接近右侧训练点后，预测带又收窄。这正是位置相关不确定性的直观表现。
这里启用了 normalize_y=True，所以核对象中的白噪声参数是在标准化后的建模尺度上优化的，learned_noise_sd=0.375 不应被直接当作原始结局尺度上的模拟噪声 0.25。predict(..., return_std=True) 返回的预测均值和标准差已经转换回原始结局尺度。
评估误差与预测不确定性
# 在测试集上获得均值和标准差，并判断观测是否落在近似95%预测区间内test_mean, test_std = gpr.predict(x_test.reshape(-1, 1), return_std=True)test_lower = test_mean - 1.96 * test_stdtest_upper = test_mean + 1.96 * test_stdis_gap = (x_test > 4.0) & (x_test < 7.0)covered = (y_test >= test_lower) & (y_test <= test_upper)# 分别汇总全体、有训练支持区和训练空白区，观察误差与不确定性是否同步增加metric_rows = []for region, mask in [(”Overall”, np.ones(x_test.size, dtype=bool)),(”Supported”, ~is_gap),(”Training gap”, is_gap),]:metric_rows.append({”Region”: region,”N”: int(mask.sum()),”RMSE”: np.sqrt(mean_squared_error(y_test[mask], test_mean[mask])),”MAE”: mean_absolute_error(y_test[mask], test_mean[mask]),”R2”: r2_score(y_test[mask], test_mean[mask]),”Coverage”: covered[mask].mean(),”MeanStd”: test_std[mask].mean(),})metrics = pd.DataFrame(metric_rows)print(metrics.round(3).to_string(index=False))# 第二张图一边检查45度一致性，一边检查预测标准差能否提示更大的绝对误差fig, axes = plt.subplots(1, 2, figsize=(9.4, 4.8))colors = np.where(is_gap, ”#D18B2C”, ”#28666E”)axes[0].scatter(y_test, test_mean, c=colors, alpha=0.56, s=24)limits = [min(y_test.min(), test_mean.min()), max(y_test.max(), test_mean.max())]axes[0].plot(limits, limits, color=”#B64D46”, linestyle=”--”, linewidth=1.4)axes[0].scatter([], [], color=”#28666E”, alpha=0.7, label=”Supported region”)axes[0].scatter([], [], color=”#D18B2C”, alpha=0.7, label=”Training gap”)axes[0].set(xlabel=”Observed outcome”, ylabel=”Predicted mean”,title=”Out-of-sample agreement”)axes[0].legend(frameon=False, loc=”upper left”)axes[0].grid(alpha=0.20)absolute_error = np.abs(y_test - test_mean)axes[1].scatter(test_std, absolute_error, c=colors, alpha=0.56, s=24)axes[1].set(xlabel=”Predictive standard deviation”, ylabel=”Absolute error”,title=”Uncertainty versus realized error”)axes[1].grid(alpha=0.20)fig.suptitle(”Gaussian process test diagnostics”, fontsize=16, x=0.06, ha=”left”)fig.tight_layout(rect=(0, 0, 1, 0.94))fig.savefig(”gpr_test_diagnostics.png”, dpi=180, bbox_inches=”tight”, facecolor=”white”)plt.close(fig)
输出：
Region N RMSE MAE R2 Coverage MeanStdOverall 240 0.259 0.206 0.909 0.962 0.347Supported 163 0.256 0.202 0.881 0.951 0.281Training gap 77 0.265 0.213 0.937 0.987 0.486
结果解释
整体测试集 RMSE 为 0.259、MAE 为 0.206、R² 为 0.909，说明后验均值较好地恢复了本次模拟中的非线性关系。左侧诊断图中，多数点靠近 45 度参考线；橙色点表示训练空白区，蓝绿色点表示有训练数据支持的区域。
更值得关注的是不确定性：有数据支持区的平均预测标准差为 0.281，而训练空白区升至 0.486，增加约 73%。空白区 RMSE 仅从 0.256 小幅升到 0.265，但区间覆盖率从 95.1% 升到 98.7%，说明模型通过扩大预测区间，对缺少训练支持作出了更保守的表达。
需要注意，98.7% 并不代表模型“比 95% 更正确”。近似 95% 预测区间依赖核函数、噪声形式和高斯假设；有限样本中的覆盖率可能高于或低于 95%。右侧图也提醒我们：预测标准差描述模型条件下的不确定性，它与单个样本最终发生的绝对误差有关，但二者不是一一对应关系。
结论
高斯过程回归的优势，不只是拟合一条灵活曲线，而是把预测值和预测不确定性放进同一套模型。只要在 predict() 中设置 return_std=True，就能进一步检查哪些区域的结论主要由数据支撑，哪些区域更多依赖核函数假设。
实际应用时还应注意三点：先对特征量纲进行处理并结合领域知识设置核；用独立测试集同时评估点预测和区间覆盖；标准高斯过程会使用全部训练样本，样本量和维度增大后计算代价会迅速上升。对于大规模数据，可考虑稀疏近似或其他可扩展的概率模型。
