高斯过程回归：小样本时间序列区间预测实战指南

做时间序列预测的人，迟早会遇到一个问题：光给一个预测值是不够的。你预测下个月销量是1000件，老板追问一句“误差多少”，你就得解释半天。更现实的情况是，供应链要备货、风控要设阈值、金融要算风险敞口，这些场景需要的不是一个孤零零的点，而是一个区间——有上界、有下界、有置信水平。这就是区间预测的价值。而在众多区间预测方法里，高斯过程回归（Gaussian Process Regression，GPR）几乎是天然为这件事设计的：它本身就能输出预测分布，也就是预测点的均值和方差，方差直接转成区间，不需要额外叠加复杂的残差分析或分位数回归。

我第一次接触GPR是在处理一批工业设备传感器的振动数据，样本量不到五十个点，跑LSTM和XGBoost都过拟合得厉害。换成高斯过程回归之后，不仅预测曲线贴合，连置信区间都有模有样，那一刻我就知道这个模型值得花时间吃透。这篇文章我会从头拆解高斯过程回归做时间序列区间预测的完整思路：核心原理、核函数选择、实操代码、预测区间构造方法、以及我在实际操作中踩过的一堆坑。适合刚接触区间预测、或者已经会用LSTM但想找小样本替代方案的读者，也适合想把预测结果从点值升级成带不确定性表达的同学。

1. 为什么区间预测比点预测更接近真实需求

1.1 点预测的局限：预测值只是条件期望

大多数时间序列模型，从ARIMA到LSTM，默认输出都是一个标量预测值。这个值的含义，本质上是给定历史数据后，目标变量条件分布的期望值。但你想想，一个分布的期望值本身携带了多少信息？它没有告诉你分布的宽度，没有告诉你尾部风险，更没有告诉你预测的可信度。比如某商品日销量过去一周是90、110、95、105、98、102、100，你预测明天是100，这个结果看起来合理。但如果过去一周的序列是50、150、40、160、45、155、100，期望同样是100，你还能淡定吗？显然不能。这两个序列的波动性完全不同，对应的预测风险也完全不一样。

点预测的问题就在这里：它把不确定性全部压缩进一个数字里，信息严重失真。在库存管理场景，点预测会直接导致安全库存设置失误；在异常检测场景，固定的点预测值会让阈值设置变成玄学；在量化交易场景，没有波动率估计的策略根本没法做仓位控制。所以很多实际项目里，真正需要交付的不是“预测值是多少”，而是“预测值大概落在哪个范围内，有多大的把握”。

1.2 不确定性来源：模型不确定性加噪声不确定性

做区间预测之前，得先搞清楚不确定性到底来自哪里。我习惯把不确定性分成两类。第一类是模型不确定性（epistemic uncertainty），指的是模型参数本身没被数据充分约束，数据不足、特征分布偏移都会放大这种不确定性。同一套数据，用不同的随机种子训练出来的LSTM权重可能差异很大，这就是模型不确定性的体现。第二类是固有噪声（aleatoric uncertainty），也就是数据本身包含的随机波动，比如传感器噪声、人为记录误差、突发因素的扰动。这类不确定性即便数据无限多也无法消除。

不同的模型处理这两类不确定性的方式差别很大。LSTM或神经网络做区间预测，通常要借助MC Dropout、Deep Ensemble或者专门的PINN损失函数来近似不确定性；分位数回归则是直接预测条件分布的分位数，本质上把不确定性内化在多个输出头里。而高斯过程回归走的是另一条路：它通过贝叶斯推断天然地同时建模这两类不确定性——模型不确定性体现在后验方差随数据稀疏程度变化，噪声不确定性则由显式的白噪声核项吸收。不需要额外搭建什么机制，预测分布直接给出完整答案。

1.3 GPR为什么天生适合区间预测

一句话概括：GPR的输出不是标量，而是分布。具体来说，对于一个测试输入$\mathbf{x}*$，GPR给出的是预测均值$\mu$和预测方差$\sigma_^2$，并且假设预测分布服从正态分布$N(\mu_, \sigma_^2)$。有了均值和方差，置信区间直接套公式：95%置信区间就是$\mu_* \pm 1.96 \cdot \sigma_*$，90%区间就是把1.96换成1.645，清清楚楚。

这和传统时间序列模型有本质区别。ARIMA也能给出预测区间，但它的区间来自残差方差的估计，隐含假设是残差独立同分布，对异方差数据无能为力。GPR的方差是输入相关的，它会在数据密集的区域给出比较窄的区间，在数据稀疏或者远离训练样本的区域自动拉宽区间。这一个性质在实际业务里极其宝贵：预测结果会自动告诉你哪些地方可信、哪些地方纯属外推，这种自适应不确定性表达，是GPR做区间预测的核心优势。

2. 高斯过程回归的核心原理拆解

2.1 高斯过程到底是一种什么东西

刚接触GPR的人很容易被“过程”这个词吓到，我换个说法你就懂了：高斯过程是对函数的概率分布。普通回归是找一个函数$f(x)$去拟合数据，GPR则是把$f(x)$本身当成一个随机变量，并且假设任意有限个输入点对应的函数值$f(x_1), f(x_2), ..., f(x_n)$都服从一个多元高斯分布。

这个假设看起来抽象，但它带来一个直接的好处：预测时不需要像神经网络那样通过梯度下降去优化一堆参数，而是可以直接用高斯分布的条件概率公式计算。给定训练数据后，任意新输入点的函数值分布可以通过协方差矩阵直接解出来。这就是为什么GPR在小样本场景非常稳，因为它的推断过程是解析的，不需要随机初始化、不需要调学习率、不需要担心梯度消失，只要核函数设定合理，结果几乎可复现。

2.2 核函数决定一切

高斯过程里最核心的组件是核函数（Kernel），或者叫协方差函数，它定义了两个样本点之间的相似度如何影响预测。核函数的本质是领域知识的注入：你觉得这个时间序列有什么性质，就用对应的核函数去表达它。

我列举几个在时间序列场景最常用到的核函数：

这里要特别强调白色噪声核的重要性。不少人第一次跑GPR，核函数只设置一个RBF，训练完发现预测曲线在训练点附近完美拟合，但预测区间窄到不真实，而且协方差矩阵经常报数值错误。原因就是没加WhiteKernel——数据里的噪声没有被显式建模，模型只能把噪声当成信号的一部分强行拟合，结果就是过拟合加数值不稳定。你把WhiteKernel加上去，噪声项会自动吸收数据中的随机波动，RBF只用去刻画真实信号的变化规律，整个模型立刻就稳了。

2.3 预测均值和方差到底怎么算出来的

这里我不堆公式推导，但建议至少理解计算流程，因为很多调参决策都建立在这个理解上。训练时，我们有训练输入$X$和观测值$\mathbf{y}$，假设潜函数$f$的先验是高斯过程，观测模型是$y = f(x) + \epsilon$，其中$\epsilon \sim N(0, \sigma_n^2)$。

训练阶段的目标是学习核函数的超参数，比如RBF的长度尺度$l$、幅度$\sigma_f^2$、噪声水平$\sigma_n^2$。怎么学？最大化边际似然（log marginal likelihood）。这个东西的作用，简单说就是评估“在当前超参数下，我观测到的数据有多大概率被这个模型生成”。通过梯度优化，可以自动找到比较合适的超参数组合。这也是GPR少有的需要迭代优化的环节。

测试阶段，给你新输入$X_$，预测分布服从$N(\mu_, \Sigma_*)$，其中：

• 均值$\mu_* = K_*^T (K + \sigma_n^2 I)^{-1} \mathbf{y}$
• 方差$\sigma_^2 = k(X_, X_) - K_^T (K + \sigma_n^2 I)^{-1} K_*$

看这两个式子你就明白，GPR预测的核心就是算一组核矩阵，然后做一次矩阵求逆。时间复杂度$O(n^3)$，这是GPR在小样本场景表现极好但在大数据集上跑不动的核心原因。区间预测需要的方差就在第二个式子——预测方差不仅取决于核函数本身，还取决于新输入和训练数据之间的协方差$K_*$。测试点离训练数据越远，核函数值越小，方差就越大。这就是自动感知数据稀疏性的机制。

3. 用GPR做时间序列预测的整体设计

3.1 数据重构：把时间序列转成回归问题

GPR本质上是个回归模型，输入输出都是向量，直接塞原始时间序列是行不通的。时间序列预测的标准做法是把序列转成监督学习格式：用过去$p$个时刻的值去预测未来$h$步的值。

假设原始序列是$x_1, x_2, ..., x_T$，要预测第$t$时刻的值，输入向量就是$[x_{t-p}, x_{t-p+1}, ..., x_{t-1}]$，输出就是$x_t$。这种嵌入方式叫滑窗法，或者叫滞后特征。窗口长度$p$的选择很重要，选太小模型学不到长周期依赖，选太大则训练样本数骤减，因为有效样本数是$T - p + 1$个。在小样本场景，这个trade-off尤其明显。我的经验是先用自相关函数(ACF)和偏自相关函数(PACF)快速看一下序列存在多少个显著的滞后期，再结合试错把$p$定在折中位置。

3.2 多步预测策略：直接法与递归法

如果你想预测未来多个时刻，不只是下一步，策略上有两种主流做法。直接法就是训练$h$个独立的GPR模型，每一个专门预测第$k$步；递归法是用一个模型，把预测值当作下一步的输入，逐步迭代出整条预测线。直接法误差不会累积，但训练成本随预测步数线性增长，而且各步之间的相关性丢失了。递归法实现简单，但训练时用的是真实历史值，预测时取而代之的是上一时刻的预测值，这个分布偏移会导致误差累积，越往后越飘。

我自己做小样本预测时更倾向递归法，因为GPR自带不确定性的特性可以弥补误差累积问题——每一步预测都把方差的传递也算进去。注意，标准的递归法其实没有显式地把方差传递到下一步输入里，严格一点的做法是采样法，即每一步从预测分布里采样一个值当作下一步输入，重复多次得到整体区间。不过这会增加计算量，小样本场景实操中直接用均值去递归也够用。

3.3 小样本场景选择GPR的深层逻辑

LSTM在小样本上的表现我已经反复提过，这里再说透一点。深度学习模型是靠大量数据驱动的，样本少时模型容量成了负担，一个几百维的LSTM隐层在几十个样本上几乎必然过拟合。正则化手段如Dropout、Weight Decay能缓解，但治标不治本。XGBoost这类树模型虽然在小样本上比深度学习稳，但它输出的依然是点预测，要搞区间预测还得额外训练分位数模型或者用残差启发式方法，两套系统的维护成本高。

GPR的优势在于它的模型复杂度可以靠贝叶斯框架自动控制。核函数里的长度尺度$l$本身就起正则化作用：数据不足以支撑复杂函数时，优化算法会把$l$调大，让函数变平滑，避免过拟合；数据充足时$l$可以变小，刻画更多细节。这种自适应正则化是GPR在小样本场景表现好的根本原因。

4. 实操流程：从模拟数据到完整预测区间

4.1 数据准备与预处理

按我的习惯，先用一个带趋势、季节性和噪声的模拟序列来演示整个流程，这样你可以直接照抄代码，再迁移到自己的数据上。

import numpy as np import matplotlib.pyplot as plt from sklearn.gaussian_process import GaussianProcessRegressor from sklearn.gaussian_process.kernels import RBF, WhiteKernel, ConstantKernel as C, ExpSineSquared from sklearn.preprocessing import StandardScaler np.random.seed(42) # 构造时间序列：趋势 + 季节 + 噪声 t = np.arange(0, 300) trend = 0.05 * t season = 5 * np.sin(2 * np.pi * t / 40) noise = np.random.normal(0, 0.8, size=len(t)) y = trend + season + noise # 训练/测试切分：用前250个点训练，预测后50个点 train_len = 250 y_train_raw, y_test_raw = y[:train_len], y[train_len:] # 标准化：GPR对量纲敏感，务必标准化 scaler_y = StandardScaler() y_train = scaler_y.fit_transform(y_train_raw.reshape(-1, 1)).ravel() y_test = scaler_y.transform(y_test_raw.reshape(-1, 1)).ravel()

标准化这一步不要偷懒。GPR的核函数超参数优化依赖对尺度比较敏感的距离计算，如果输入输出量纲差距很大，要么初始核参数选不好，要么优化陷入糟糕的局部最优。用StandardScaler把数据归一化到零均值单位方差，是整个流程里性价比最高的一步。

4.2 构造滑窗数据集与核函数初始化

def create_sequences(y_data, window_size): X, Y = [], [] for i in range(window_size, len(y_data)): X.append(y_data[i-window_size:i]) Y.append(y_data[i]) return np.array(X), np.array(Y) window = 20 X_train, Y_train = create_sequences(y_train, window) X_test_raw_seq, Y_test_true = create_sequences(np.concatenate([y_train, y_test]), window) X_test = X_test_raw_seq[train_len-window:] # 对齐测试集起点 Y_test_true = Y_test_true[train_len-window:]

窗口大小取20，也就是用过去20个时刻预测当前时刻。这个选择兼顾了趋势和40步周期的季节性模式，20个点足够覆盖半个季节周期。如果数据周期更长，这里也要相应拉大。

然后是核函数配置。以这个模拟序列为例，我选用组合核：

kernel = C(1.0, (1e-4, 1e2)) * RBF(length_scale=10.0, length_scale_bounds=(1e-2, 1e3)) + C(1.0, (1e-4, 1e2)) * ExpSineSquared(length_scale=10.0, periodicity=40.0, periodicity_bounds=(35, 45)) + WhiteKernel(noise_level=0.1, noise_level_bounds=(1e-5, 1e1)) gpr = GaussianProcessRegressor(kernel=kernel, alpha=1e-8, normalize_y=True, n_restarts_optimizer=5) gpr.fit(X_train, Y_train)

这套核函数配置我拆开讲一下。第一个部分是常数核乘以RBF，用来刻画整体的平滑趋势和局部变化；第二部分是常数核乘以周期核，用来刻画40步的季节周期，我显式限制了periodicity在35到45之间，防止优化器把周期调到离谱的值；最后是白噪声核，吸收观测噪声。很多初学者只用一个RBF，结果预测结果变成一条插值曲线，问题就出在缺少结构化的先验。时间序列数据里，趋势和季节性往往是同时存在的，核函数至少要表达这两类特征，否则再强的拟合能力也是白搭。

4.3 训练、预测与区间构造

gpr.fit(X_train, Y_train) y_pred_mean, y_pred_std = gpr.predict(X_test, return_std=True) # 反标准化 y_pred_mean = scaler_y.inverse_transform(y_pred_mean.reshape(-1, 1)).ravel() y_pred_std = y_pred_std * scaler_y.scale_.ravel()[0] # 方差随缩放等比例变换 # 95%置信区间 lower = y_pred_mean - 1.96 * y_pred_std upper = y_pred_mean + 1.96 * y_pred_std

注意反标准化的时候要把标准差也做逆变换。标准化后的预测方差是标准正态空间里的方差，乘以原始尺度的scale值就能还原到原始量纲。这里容易出错，很多人把均值和方差混在一起直接inverse_transform，结果区间要么大得离谱要么窄得可笑。

画出来看效果：均值和真实曲线基本重叠，区间在数据模式比较固定的地方收窄，在趋势出现转向或者测试集较远的位置变宽。这个区间本身就告诉你模型的自信程度，不用你去猜。

4.4 区间质量怎么量化评价

光画图还不够，项目交付需要量化指标。两个指标基本够用：预测区间覆盖率（PICP）和区间平均宽度（PINAW）。PICP衡量真实值落在预测区间内的比例，公式是：

coverage = np.mean((Y_test_true >= lower) & (Y_test_true <= upper))

PINAW则是区间宽度的平均值除以数据范围，用来评估区间是否过宽：

pinaw = np.mean(upper - lower) / (np.max(y_test_raw) - np.min(y_test_raw))

理想的预测区间是：覆盖率接近置信水平（比如95%区间就接近0.95），同时区间宽度尽可能窄。如果覆盖率远低于0.95，说明模型过度自信，区间太窄；如果覆盖率远超0.95，PINAW又很大，说明区间太保守，预测信息量低。这两个指标要放在一起看，不会有人只要覆盖率不要宽度的。

4.5 用真实数据时的替代方案

如果是小样本仿真数据或者公开数据集，直接把上面的流程照搬就行。但真实业务数据通常有缺失值和异常点。缺失值处理上，GPR本身不能天然处理NaN，滑窗构造前需要用线性插值或前向填充补全。异常点处理要谨慎，GPR的噪声核能吸收一定程度的离群值，但严重异常值会扭曲核超参数的学习，建议先用IQR或MAD等稳健统计方法标记异常点，再决定是剔除还是缩尾处理。

5. 实操中遇到的典型问题和排查方法

5.1 预测区间窄得不正常或者宽得离谱

区间过窄，先检查数据标准化是否做了、WhiteKernel是否加了、训练数据是否过少导致核函数过度自信。区间过窄最典型的原因是噪声核的初值设得太小，优化器找不到更大噪声水平的最优解，模型把噪声当信号拟合，方差被严重低估。解决办法是放宽噪声核的搜索范围，比如把noise_level_bounds的下界降到1e-5，上界提升到100，给优化器更多选择余地。

区间过宽的情况多发生在测试集远离训练集时。这是GPR的正常行为，因为外推区域的方差本就该大。但如果离训练集很近的区间也很宽，多半是核函数的长度尺度太短，使得样本点之间的相关性很弱，模型无法从邻近的信息中借力。此时适当增大RBF的length_scale初始值，或者检查是否有周期性核却设成了很短的周期。

5.2 训练报错：协方差矩阵不满足正定

这个错误我刚开始跑GPR时频繁遇到。报错信息通常是“Cholesky decomposition failed”或者“singular matrix”。根本原因是核矩阵的条件数过大，接近奇异，也就是矩阵里有几乎线性相关的行。触发情况主要有三种：训练样本里有重复值或近似重复值导致核矩阵对角线以外的元素趋近于1；噪声核的噪声水平设置得过小，对角线上的微小扰动不够打破病态；数据本身有大量相同的输入组合，比如某段时间传感器一直归零。

排查顺序也是先加/调大WhiteKernel的噪声水平，再检查训练样本是否有重复或近似重复的情况，必要时做去重处理，最后在GaussianProcessRegressor里设置alpha参数，比如1e-8，相当于给对角线加一个小的正则化项。我自己一般三者叠加使用，基本能覆盖所有场景。

5.3 超参数优化陷入局部最优

Sklearn的GPR通过优化边际似然来学习超参数，但边际似然函数在参数空间中是非凸的，很可能有多个局部极大值。直接跑默认参数几乎必然被坑。两个手段解决：一个是设置n_restarts_optimizer，这个参数让优化器从多个随机初始点出发去搜索，我一般设5到10，效果立竿见影；另一个是手动约束核函数的搜索范围，比如周期核的periodicity_bounds，如果你知道业务周期不是每12个月，就老老实实把范围限制在合理区间，别指望优化器从完全随机的初值里找到正确答案。

还有个小技巧，训练前打印一下优化后的超参数值，如果发现某些参数贴在边界上，大概率是初始值不合适或者边界设置太紧，调整边界重新训练。

5.4 时间序列场景特有问题

递归多步预测时，误差会随着步数累积。如果你用递归法往后预测50步，最后几步的区间宽到没有业务价值，这不是GPR的缺陷，而是递归策略的天性。缓解办法是改用直接法做长程预测：每步用独立训练的模型，只是工作量会翻倍。

另一个常见问题是数据的平稳性。GPR本身不对序列的平稳性做硬性要求，核函数里的趋势项可以捕捉线性趋势，但如果序列有指数级增长或者明显的结构性突变，GPR就会表现得力不从心。一种做法是先做差分，用差分后的平稳序列训练GPR，再把预测结果累加回去；另一种做法是用对数变换压缩量纲。我在业务中遇到不同股票价格序列时，通常两种方法配合使用，效果比直接用原始序列好很多。

5.5 GPR和LSTM怎么选

很多人看到“LSTM时间序列预测python”就冲进深度学习，但选模型要看你手里有什么牌。样本量在一两百以下、特征维度不高、预测需要带不确定性表达，选GPR基本是毋庸置疑的；数据量大、模式极其复杂、非线性程度高、样本上千起步，再考虑LSTM。

这里给一个粗参考线：如果你手头数据量小于200条，别犹豫，直接GPR；200到1000条试着用简化版LSTM配合强正则化，但也要准备好区间预测的后处理方案；上千条后再用深度模型。时间序列领域有个悖论——“深度学习很强大”，但大部分实际业务场景根本没有足够的历史数据去喂养它。GPR作为小样本区间预测的答案，很多时候不是备选，而是正解。

6. 进一步扩展的空间

做完了基础的区间预测之后，这套框架还可以继续往几个方向深入。一个是多输出GPR，针对多个相关时间序列可以联合建模，利用变量间的相关性生成更紧密的预测区间。另一个是异方差噪声的建模，标准GPR假设噪声水平全局一致，但真实数据的波动往往在特定时段变大，比如交易时段的波动率显著高于休息时段，这是异方差现象。处理这类问题，可以借助稀疏高斯过程或者变分推断方案，用输入相关的噪声函数替代固定白噪声核，让区间宽度随输入变化更贴合实际波动模式。

还有异常检测场景。GPR预测区间的上下界天然就是动态阈值，突破区间的点可以标记为异常。相比传统的固定阈值或者滑动窗口阈值，GPR阈值的优势在于它会基于输入的局部结构自适应调整宽度，不会在数据波动本来就大的地方频繁误报。我在传感器故障检测里试过这个方案，误报率比固定阈值下降了一个档次。

最后说一下我个人的选择偏好。我在实际项目里最常用的核函数组合是：常数核乘RBF加周期核加白噪声核，这套组合覆盖了趋势、季节性和噪声三大要素，适配绝大多数具有周期性的业务数据。如果没有明显的周期特征，就把周期核删掉，保留RBF和白噪声核。跑通流程之后再根据数据特征逐步加复杂度，不要一开始就堆一大堆核函数——GPR的核函数组合是乘法关系，核越多，优化空间就越复杂，局部最优问题会更严重。先把基本盘做稳，再谈扩展。

关于预测区间的业务落地，我还有一个屡试不爽的经验：不要把95%区间直接甩给业务方，而是同时提供均值、80%区间和95%区间三层结果。80%区间给日常决策用，95%区间给风险控制场景用，均值作为最可能的估计，这样业务方既不会觉得区间太宽没用，也不会在风险场景里因为区间过窄而误判。这套区间预测框架的意义，不只是换了个算法，而是把预测从一个数字变成了一套带有可靠性量化的决策工具。