GP-EnKF在线高斯过程回归：把核超参数塞进集合卡尔曼滤波

简介：GP-EnKF 是一份面向在线高斯过程回归学习与实现的 Python 代码资源，对应 Fusion 2018 论文思路，适合已具备高斯过程、卡尔曼滤波基础的研究生、算法工程师和科研人员，用于理解并复现 EnKF 如何在线更新高斯过程先验、完成实时预测与不确定性估计。资源以 zip 压缩包形式提供，大小约 22KB；上游未给出有效文件总数与类型明细，包体较精简，重点应落在核心实现代码与论文相关示例上，便于按模块阅读和二次改造。页面显示已有 307 人学习或下载。读者可从中获得 GP-EnKF 的预测—更新循环实现线索，包括集合状态构建、观测融合、超参数与核函数设置等关键环节，并参考环境科学、控制工程、信号处理等场景，将在线数据融合、动态系统状态估计与概率预测方法迁移到自己的实验中。

1. GP-EnKF 代码包拆解：在线高斯过程回归为什么要把超参数塞进卡尔曼滤波

传感器每 200 ms 吐一个点，模型必须在这个窗口里把核超参数和预测分布一起更新完，还不能把历史数据全留着——这是在线高斯过程回归最真实的工况。GP-EnKF 这份 Python 代码处理的就是这个场景：它把 GP 的诱导点函数值和核超参数打包成一个状态向量，用集合卡尔曼滤波做逐点的递推贝叶斯更新，思路来自 FUSION 2018 那篇在线 GP 回归与学习的滤波工作。它更像一份可以直接投入使用的在线推断基础框架，不需要每次进来新数据就重训，也不用对边缘似然求梯度；它并非 GP 库的又一层封装。

适合三类人：手里有实时回归任务、又必须要不确定度输出的工程师；给 MPC、状态估计、软测量找在线代理模型的开发者；以及读过高斯过程教材、但对「在线」始终没找到能跑起来的参考实现的同学。

代码本身不长，真正费时间的是几个参数和数值细节：集合开多大、诱导点怎么摆、观测噪声估多少、协方差要不要膨胀。后面按这条线拆。

2. 从批量 GP 到在线 GP：把 O(N³) 的账算清再谈滤波

2.1 批量 GP 卡在哪一步：不是预测慢，是每次都要重来

标准 GP 回归的后验均值是 $\mu_* = k_^\top (K+\sigma_n^2 I)^{-1} y$，预测方差是 $k_{**} - k_^\top (K+\sigma_n^2 I)^{-1} k_*$。整套计算的核心是一次 Cholesky 分解，复杂度 O(N³)。N=2000 时在普通笔记本上大概几百毫秒量级，具体数字跟 BLAS 后端和 CPU 关系很大，我建议你自己测一遍——重点不在于绝对值，而在于增长率：数据翻倍，时间大约涨 8 倍。

真正致命的是超参数。边缘似然每评估一次就要一次 Cholesky，L-BFGS 迭代几十步、每步还要算梯度（梯度里同样含 $K^{-1}$），一轮训练就是几十次 O(N³)。数据以 5 Hz 进来，一小时 18000 条，批量重训这条路径在实时系统里彻底走不通。

更麻烦的是「重训」这件事本身的语义：新数据进来重新优化超参数，解会跳。上一秒预测方差还是 0.02，下一秒变成 0.15，下游的控制器会以为模型崩了。在线 GP 要解决的，不只是计算量的问题，还有解路径的连续性。

2.2 三条主流在线化路线的取舍

滑动窗口最省事，代价是超参数在窗口切换时跳变；增量 Cholesky（把 $K+\sigma^2 I$ 的分解做 rank-1 更新）速度快，但超参数被钉死，一旦工况漂移，模型不会自己调长度尺度，只能靠人工定期重训——这就是「稀疏 GP 上线三个月后精度莫名下滑」这类问题的常见根源。

滤波路线把 GP 的诱导点参数化写成状态空间模型，用卡尔曼类方法递推。它的好处是把「学超参数」和「估后验」变成同一件事：都是状态更新。代价是数值细节多，尤其是协方差的秩和膨胀。

2.3 为什么是 EnKF，而不是 EKF 或 UKF

观测算子 $h(\cdot)$ 里含核矩阵求逆，而核矩阵又依赖状态里的长度尺度，所以 $h$ 对状态是非线性的。用 EKF 要推 $\partial h/\partial \theta$，也就是 $\partial (K_{zz}^{-1} k_*) / \partial \log \ell$，链式展开下来一长串，且每步都要算一次矩阵导数，工程量与出错概率都高。UKF 不用雅可比，但 sigma 点数量是 $2n+1$，状态维度 $n = M+2$ 时（诱导点 40 个就是 84 个 sigma 点）每步要算 84 次核矩阵分解，成本不划算。

EnKF 用随机采样的集合近似协方差，只需前向计算 $h$，天然导数-free，集合大小固定在 80~100 就能工作，跟状态维度基本解耦。这就是这份代码选 EnKF 的核心动机——工程上可维护，而且超参数能跟着一起更新。

2.4 先把批量基准跑通，再谈在线

任何在线实现都需要一个离线基准来对齐，否则你没法判断在线版的偏差是算法问题还是实现 bug。下面这段是批量 GP 的最小实现，后面所有在线结果都拿它对照。

import numpy as np from scipy.linalg import cho_factor, cho_solve def rbf(A, B, ls, sf): """各向同性 RBF 核。A:(n,d) B:(m,d) -> (n,m)""" d2 = (A * A).sum(1)[:, None] + (B * B).sum(1)[None, :] - 2.0 * A @ B.T d2 = np.maximum(d2, 0.0) # 浮点误差可能给出 -1e-16 return sf ** 2 * np.exp(-0.5 * d2 / ls ** 2) def gp_batch(X, y, Xs, ls=1.0, sf=1.0, sn=0.1): """批量 GP 后验预测，返回均值与方差""" n = len(X) K = rbf(X, X, ls, sf) + (sn ** 2 + 1e-8) * np.eye(n) # 1e-8 是 jitter c = cho_factor(K, lower=True) alpha = cho_solve(c, y) # (K+σ²I)^{-1} y Ks = rbf(Xs, X, ls, sf) # (m, n) mu = Ks @ alpha v = cho_solve(c, Ks.T) # (n, m) var = sf ** 2 - np.einsum('ij,ji->i', Ks, v) # 逐点取对角 return mu, np.clip(var, 1e-12, None) # 数值上可能微负

cho_factor返回的是下三角分解结果，cho_solve复用同一份分解同时算均值和方差——这是批量 GP 唯一值得强调的性能细节，别写成两次np.linalg.solve。jitter加在训练协方差对角上，1e-8 起步，如果你的长度尺度很小或诱导点很密，这个值要往上调到 1e-6。

拿这段代码在一段模拟数据上跑一遍、记下耗时，你就有了对照组。我一般会把 N 从 200 拉到 3200 做一条耗时曲线，确认它确实接近立方增长——这条曲线是后面判断在线版「省了多少」的唯一依据。

3. 把 GP 后验塞进 EnKF：状态向量、观测算子与集合更新

3.1 状态向量怎么定：诱导点值在前，对数超参数在后

GP 的诱导点参数化（sparse pseudo-input 那一套）把函数写成 $f(x) = k(x,Z)K_{zz}^{-1}u$，其中 $Z$ 是 $M$ 个诱导点位置，$u$ 是这些位置上的函数值，先验 $u \sim \mathcal{N}(0, K_{zz})$。这把一个无穷维的函数推断压缩成了 $M$ 维向量的推断，而 $M$ 可以直接作为状态维度。

状态向量定义为：

$$z = [,u_1, \dots, u_M,\ \log \ell,\ \log \sigma_f,]^\top$$

长度 $n_{state} = M + 2$。用对数参数化而不是原始值，是因为长度尺度和幅值必须为正，而 EnKF 的更新是线性叠加、不保号，直接更新 $\ell$ 很容易跑出负值然后核矩阵集体失效。取对数之后，状态可以自由取任意实数，指数变换回去恒正。这是整个实现里最便宜也最关键的一个设计。

$u$ 的先验直接从 $\mathcal{N}(0, K_{zz})$ 采样，集合初始化时每个成员独立抽一次。超参数则围绕一个初值加小扰动，扰动幅度取 0.05~0.1（对数尺度），太大集合一开始就散掉，太小则初始协方差接近奇异、第一次更新几乎不动。

提示：诱导点位置 $Z$ 不放进状态向量。它们是「支撑点」，不是待估量；一旦允许它们动，状态维度会随 $M \cdot d$ 暴涨，且 $Z$ 变化会让 $K_{zz}$ 整体改变，集合的一致性没法维持。实践里 $Z$ 固定，只在离线阶段选。

3.2 观测算子与扰动观测：EnKF 里最容易被写错的一行

观测算子 $h(z, x)$ 把状态映射到「在输入 $x$ 处的预测值」：

import numpy as np from scipy.linalg import cho_factor, cho_solve def h_member(u, log_ls, log_sf, xq, Zind, jitter=1e-6): """单个集合成员的观测预测：k(x*,Z) K_zz^{-1} u""" ls, sf = np.exp(log_ls), np.exp(log_sf) Kzz = rbf(Zind, Zind, ls, sf) + jitter * np.eye(len(Zind)) c = cho_factor(Kzz, lower=True) kx = rbf(xq[None, :], Zind, ls, sf).ravel() # (M,) w = cho_solve(c, kx) # K_zz^{-1} k(x*,Z) return float(w @ u) def H_ensemble(Zf, xq, Zind): """把整个集合映射成观测预测集合。Zf:(n_state, n_ens) -> (n_ens,)""" M_ind = len(Zind) return np.array([ h_member(Zf[:M_ind, i], Zf[M_ind, i], Zf[M_ind + 1, i], xq, Zind) for i in range(Zf.shape[1]) ])

注意h_member里每个成员用的是自己那份超参数去算 $K_{zz}$，不是共享一份。这一点写错的话，超参数就退化成「只影响输出尺度」的装饰，学不出任何东西——因为集合里所有成员共用同一个核矩阵，超参数只在最后乘 $u$ 时才起作用，梯度信息全丢。

分析步按标准 EnKF 写：

def enkf_analysis(Zf, y, Hf, R, rng, inflation=1.02): """ Zf : (n_state, n_ens) 预报集合 y : float 当前观测 Hf : (n_ens,) 观测预测集合 R : float 观测噪声方差 """ n_ens = Zf.shape[1] # 协方差膨胀：先把集合往均值外侧推，抵消有限的秩导致的信息丢失 mean = Zf.mean(axis=1, keepdims=True) Zf = mean + inflation * (Zf - mean) Hf = Hf.mean() + inflation * (Hf - Hf.mean()) dz = Zf - Zf.mean(axis=1, keepdims=True) # (n_state, n_ens) dh = Hf - Hf.mean() # (n_ens,) Cyy = dh @ dh / (n_ens - 1) + R # 标量 Cxy = dz @ dh / (n_ens - 1) # (n_state,) K = Cxy / Cyy # 增益向量 # 每个成员一份独立的观测扰动，方差必须等于 R D = y + rng.normal(0.0, np.sqrt(R), size=n_ens) Za = Zf + np.outer(K, D - Hf) return Za, K

D那一行是整个 EnKF 最容易出问题的地方。观测扰动必须逐成员独立、方差严格等于 $R$。如果所有成员共用同一个扰动，观测噪声的影响会被系统性抵消，集合方差被低估，跑几百步之后集合直接坍缩成一条线；如果扰动方差写成 $R/2$ 或者忘了加，模型会变得过度自信，预测区间窄得离谱。

协方差膨胀inflation默认 1.02，是在补偿「集合大小有限 → 样本协方差秩只有 n_ens-1 → 低估真实不确定性」这个结构性偏差。数据集短、观测密集时可以开到 1.05；开到 1.1 以上通常表现为预测方差虚高、均值抖动。

3.3 分析步与预报步：一个循环里的两次更新

EnKF 的一步分两段。分析步吃进当前观测，把集合往观测拉；预报步把状态推进到下一个时刻。

预报步里，$u$ 是诱导点上的函数值，属于模型的静态参数，理论上不随时间变，所以用持续性预报 $u_{k+1} = u_k$。超参数不一样：工况在漂，长度尺度应该缓慢跟着走，用随机游走：

$$\log \theta_{k+1} = \log \theta_k + \eta, \qquad \eta \sim \mathcal{N}(0, q^2)$$

def forecast_step(Za, rng, rw=1e-3): """持续性预报：u 不动，对数超参数加随机游走""" Zf = Za.copy() n_state = Zf.shape[0] Zf[n_state - 2:, :] += rng.normal(0.0, rw, size=(2, Zf.shape[1])) return Zf

rw这一项是整个实现里唯一调节「超参数适应速度」的旋钮，量纲是每次更新的对数标准差。1e-3 意味着每步长度尺度变化约 0.1%，1000 步累积约 3 倍——对小时级的缓慢漂移是合适的。数据是平稳的就把rw压到 1e-4 甚至 0，让超参数几乎不动；工况切换频繁（比如设备启停）时抬到 5e-3，但代价是超参数会一直小幅抖动，预测方差偏高。

注意这里给超参数加了过程噪声，却没给 $u$ 加。有些实现为了维持集合离散度也会给 $u$ 加一点，我不建议——$u$ 的先验方差本来就由 $K_{zz}$ 决定，人为加噪声等于篡改先验，会让预测方差的校准彻底失真。

3.4 和标准卡尔曼更新式的对照

把两套公式并排放，能看清 EnKF 到底「近似」掉了什么：

关键差別在最后一行：EnKF 不维护 $P$，它维护集合本身，协方差是从集合现算的。这带来两个直接后果——第一，$K$ 的秩最多 $n_{ens}-1$，状态维度超过这个数时增益必然是降秩的；第二，更新后的协方差自动包含非线性传播的效果，不需要雅可比。所以 $n_{ens}$ 一般取 80~100，而 $n_{state}$ 最好控制在 60 以内，超了就上 localization 或者减诱导点。

4. 诱导点的在线维护：代码骨架、初始化与四个必调参数

4.1 诱导点从哪来：三种摆法及各自的失效场景

诱导点的位置决定了模型能表达什么。常见的三种选法：

均匀网格。1D 或 2D 输入下最省事，直接按分位数切。缺点是数据分布不均匀时，密集区采样不足、稀疏区浪费预算。我一般先用它跑通流程，再换。

k-means 聚类中心。对训练输入跑一次 k-means，取 $M$ 个中心。这是覆盖数据密度最稳妥的做法，代价是每来一批新数据要重跑聚类——但只在离线阶段做一次，在线阶段固定，成本可以接受。注意 k-means 对量纲敏感，先做标准化。

在线增删。新数据点离最近诱导点的距离超过阈值 $\tau$ 时把它加进去，预算超了就合并距离最近的一对。听起来很美，实际有两个硬伤：状态向量长度一变，整个集合必须重新初始化，前面学到的后验全丢；合并操作让 $u$ 的语义发生跳变。绝大多数能稳定跑起来的实现都固定 $Z$，把「覆盖不足」的问题交给长度尺度和更大的 $M$ 来吸收。

下面这个初始化函数用的是 k-means 路线，没有 sklearn 依赖时可以用分位数采样替代：

def init_inducing(X, M, rng, mode="quantile"): """从训练输入里挑 M 个诱导点位置""" if mode == "quantile": # 每个维度独立取分位数，再组合。1D/2D 效果好，高维退化为网格 qs = np.linspace(0.05, 0.95, M) cols = [np.quantile(X[:, j], qs) for j in range(X.shape[1])] Z = np.stack(cols, axis=1) if X.shape[1] == 1: Z = Z.reshape(-1, 1) return Z # 随机抽样版本：从数据点里无放回抽 idx = rng.choice(len(X), size=M, replace=False) return X[idx].copy()

分位数版本在 1D 下等价于等间距网格，在 2D 下会退化成对角线上的一条线——这是它的硬伤，2D 以上请老实上 k-means。高维输入（$d>5$）时诱导点数量要按 $M \approx 10d$ 起步，否则 $K_{zz}$ 条件数会很差。

4.2 预测均值与方差的集合估计

给定训练好的集合，在查询点 $x_*$ 处的预测就是逐成员算一遍再统计：

def predict(Za, xq, Zind, sn=0.1, ddof=1): """集合预测：返回均值、潜函数方差、含噪方差""" M_ind = len(Zind) n_ens = Za.shape[1] vals = np.array([ h_member(Za[:M_ind, i], Za[M_ind, i], Za[M_ind + 1, i], xq, Zind) for i in range(n_ens) ]) mu = vals.mean() var_f = vals.var(ddof=ddof) # 潜函数方差 return mu, var_f, var_f + sn ** 2 # 后者是含观测噪声的预测方差

ddof=1不能省。EnKF 里的样本协方差全部用无偏估计，混用ddof=0会让方差系统性偏小 $1/n_{ens}$，跑长序列时这个偏差会累积成明显的过度自信。

潜函数方差和含噪方差要分开返回。下游如果拿预测做采集函数或者置信区间判断，用的是潜函数方差；如果做的是传感器读数预测，用含噪方差。这两个混用是最常见的接口层 bug——表现是预测区间看起来「刚好」，但覆盖率测试只有 70%。

4.3 在线主循环与参数表

把上面几块拼起来就是完整的在线循环。注意初始化顺序：先采 $u$，再采超参数，两者独立。

def run_online(X, y, Zind, n_ens=100, sn=0.1, rw=1e-3, inflation=1.02, seed=0, warmup=30): rng = np.random.default_rng(seed) M_ind = len(Zind) n_state = M_ind + 2 # ---- 集合初始化 ---- Kzz0 = rbf(Zind, Zind, 1.0, 1.0) + 1e-6 * np.eye(M_ind) L = np.linalg.cholesky(Kzz0) U = L @ rng.normal(size=(M_ind, n_ens)) # u ~ N(0, K_zz) log_ls = np.log(1.0) + rng.normal(0, 0.05, size=n_ens) log_sf = np.log(1.0) + rng.normal(0, 0.05, size=n_ens) Z = np.vstack([U, log_ls[None, :], log_sf[None, :]]) hist = [] for k in range(len(X)): # 预热期只用先验，不做更新，让集合先散开 if k >= warmup: Hf = H_ensemble(Z, X[k], Zind) Z, K = enkf_analysis(Z, y[k], Hf, sn ** 2, rng, inflation) hist.append(Z.copy()) Z = forecast_step(Z, rng, rw) return Z, hist

warmup这一段容易被忽略。前几十步如果不做更新、只走预报，集合会按先验方差自然散开；如果从第一步就开始更新，初始集合的离散度被人为设得很小（0.05 的对数扰动），第一次更新几乎不动，之后要很久才能恢复正常的离散度。我一般设warmup=30，数据量少于 100 条时干脆关掉在线更新，直接当批量 GP 用。

四个必须调的参数，连同调错的症状：

sn最容易被随手填成 0.1。它的正确量级应该来自传感器手册或者一段静态数据的残差标准差。填大了模型反应迟钝，填小了每次更新都当成真信号去追，超参数会被噪声推得到处跑。

5. 踩坑排查：集合坍缩、负方差与超参漂移的五个现场

5.1 集合坍缩：预测方差一路掉到 1e-8

现象：跑了两三百步之后，predict返回的var_f单调下降，最后稳定在 1e-8 量级，预测区间变成一条线，但 RMSE 并没有同步变好。

原因：EnKF 的更新是收缩性的，每次分析步都会把集合往观测方向压，样本协方差单调减小。有限集合（$n_{ens}=100$）下这个衰减比理论值快得多，因为增益的秩只有 99，一部分方向上的方差被永久抹掉。没有膨胀机制的话，几百步之后就彻底坍缩。

解决：先确认inflation真的生效了——把它放在enkf_analysis内部、算协方差之前，而不是在调用处乘。然后观察集合标准差的时间序列，正常应该稳定在一个平台上小幅波动。如果 1.02 压不住，往上试 1.05；还压不住，说明n_ens相对n_state太小，加集合或者减诱导点。我遇到过一次 $n_{state}=82$ 配 $n_{ens}=50$ 的组合，膨胀开到 1.08 才勉强稳住，最后还是把诱导点砍到 30 才正常。

5.2 负方差与 Cholesky 报错：LinAlgError: Matrix is not positive definite

现象：某一步之后cho_factor抛异常，或者预测方差出现 -1e-9 这样的负值被clip掉了。

原因：两个来源。一是Kzz + jitter*I在超参数被更新到极小长度尺度时条件数爆炸，1e-6 的 jitter 不够；二是集合里出现了几乎重复的成员（比如某些成员的 $u$ 被更新到同一值），样本协方差奇异。

解决：jitter 改成自适应——先算np.linalg.cond(Kzz)，超过 1e10 就把 jitter 从 1e-6 抬到 1e-4；同时给 $\log \ell$ 加硬下界，更新后np.clip(log_ls, np.log(0.01), np.log(100))。这个 clip 看着粗暴，但它挡住的是「长度尺度跑到 1e-5 导致核矩阵退化成单位阵」这类灾难，比事后调试便宜得多。负方差用np.clip(var, 0, None)兜底，同时打日志——出现频率高说明模型已经不对了，别当成正常现象。

5.3 超参数漂移：长度尺度单调爬到边界

现象：$\log \ell$ 的时间序列单调上升或下降，最后贴着 clip 的边界不动，模型退化成常数预测或者纯噪声拟合。

原因：rw太大加上sn估得太小。观测噪声小 → 每次更新的增益大 → 超参数被单个观测点猛推一把 → 随机游走又把它往同方向再推一点，形成正反馈。

解决：先把sn按实际噪声标定，再把rw砍一个数量级观察。另一个有效的做法是给超参数的随机游走加一个弱回归项（往初始值拉），相当于给超参数一个高斯先验：

def forecast_step_anchored(Za, rng, rw=1e-3, anchor=0.0, kappa=1e-3): Zf = Za.copy() n = Zf.shape[0] Zf[n-2:, :] += rng.normal(0.0, rw, size=(2, Zf.shape[1])) Zf[n-2:, :] -= kappa * (Zf[n-2:, :] - anchor) # 弱回归，抑制漂移 return Zf

kappa取rw的十分之一左右。它不改变稳态性能，但能把发散的时间尺度从几百步拉到几千步。

5.4 诱导点数超过集合数：增益矩阵被悄悄降秩

现象：模型能跑，不报错，但拟合效果始终比同配置的批量 GP 差一截，残差里能看到明显的系统性结构（不是白噪声）。

原因：$n_{ens} < n_{state}$ 时，样本协方差矩阵 $\hat P$ 的秩最多 $n_{ens}-1$，增益 $K$ 把状态投影到了一个低维子空间。落在零空间里的方向永远得不到更新，等价于那部分诱导点的函数值被冻结在先验上。$M=100$、$n_{ens}=100$ 就是典型的踩线配置。

解决：保持 $n_{ens} \ge 2 n_{state}$ 是安全的，做不到就减 $M$。另一个办法是给协方差做局部化（localization），但 GP 的状态没有空间索引，做不了标准的距离衰减局部化；实践中更简单的是把诱导点分成几组、轮流更新，或者干脆接受更大的集合。跑之前先打印一下n_state和n_ens，这两个数看一眼就能避免大部分问题。

5.5 数据顺序敏感：同一批数据换个顺序结果差很多

现象：把同一段数据打乱顺序喂进去，最终的 RMSE 差 20% 以上，超参数收敛到完全不同的值。

原因：在线递推本身就是路径依赖的——EnKF 只保留当前集合，不保留历史数据，所以信息是有损的。顺序影响来自两方面：观测噪声被假设为各向同性同方差，而实际噪声可能有时变特性；超参数的随机游走让早期数据的影响被逐渐遗忘（这正是在线学习的本意，但遗忘速率由rw决定）。

解决：先接受它。路径依赖是在线 GP 的固有性质，不是 bug。要减小它，就降rw、加warmup、把inflation调小。要量化它，就跑 3~5 个不同的打乱顺序，看指标的标准差——如果标准差超过均值的 10%，说明rw太大或者数据量对当前M来说不够。我一般会在评估报告里同时给「顺序 A」和「顺序 B」两组数，而不是只报最好那一组。

6. 收敛验证：三条曲线判断在线 GP 是否真的在学

在线模型最麻烦的地方是它永远不「训练完」，你没法用最终 loss 判断好坏。我的做法是固定看三条曲线，缺一条就不敢上线。

第一条：一步预测的标准化残差。定义 $e_k = (y_k - \mu_k) / \sqrt{\sigma_k^2 + \sigma_n^2}$，正常应该近似标准正态——均值在 0 附近、标准差在 1 附近、无自相关。如果标准差到了 3，说明模型过度自信（sn太小或者inflation不够）；如果只有 0.3，说明区间太宽（sn太大或者rw太大）。这条曲线比 RMSE 敏感得多，RMSE 对区间宽度完全不敏感，而在线 GP 的价值恰恰在区间上。

第二条：集合离散度与残差的一致性（spread-skill）。把窗口内所有时刻的集合标准差求平均，再除以同期残差的 RMS，比值应该落在 0.8~1.2。低于 0.5 就是坍缩的前兆，高于 2 说明集合被人为撑得太大。这个指标能提前几百步预警 5.1 里的坍缩问题。

第三条：超参数轨迹。把 $\log \ell$ 和 $\log \sigma_f$ 的时间序列画出来。健康的轨迹是「快速调整一段 → 趋于平稳 → 小幅波动」，如果一直在单调爬升，说明 5.3 的漂移已经发生；如果从第一步到最后一步几乎不动，说明rw压得太死，超参数实际上没在学，退化成固定超参的稀疏 GP。

有一个小技巧值得单独说：用集合本身做交叉验证。集合里每个成员都可以看作后验的一个样本，把第 $i$ 个成员单独拿出来算它在前 $k$ 步的预测 NLL，取平均，就得到一个不需要额外留出集的泛化估计。代价是 $n_{ens}$ 次前向计算，但比重新跑一遍在线流程便宜得多。我在数据量紧张的场景下经常用它来做超参数的粗筛：

def loo_nll(Za, X, y, Zind, sn=0.1): """用集合成员近似后验样本，算平均 NLL""" n_ens = Za.shape[1] nlls = [] for i in range(n_ens): col = Za[:, i:i+1] # 单个成员当作一个后验样本 err = 0.0 for k in range(len(X)): mu, var_f, var_y = predict(col, X[k], Zind, sn) err += 0.5 * np.log(2 * np.pi * var_y) + 0.5 * (y[k] - mu) ** 2 / var_y nlls.append(err / len(X)) return float(np.mean(nlls)), float(np.std(nlls))

返回的标准差本身也有信息量：它反映后验的分散程度，标准差大说明不同成员给出的解释差很多，模型还没收敛。

从那以后我每次上线一个新的在线 GP 模型，都强制先跑一遍这三条曲线加一次 LOO-NLL，哪怕指标看着漂亮也不跳。在线模型的失败几乎都是缓慢的、不报错的——它不会抛异常，只会安静地变得过度自信，然后在某个工况切换的瞬间给出一个窄得离谱的预测区间。这三条曲线就是那个后悔药。希望帮到你。

本文还有配套的精品资源，点击获取