CUSUM 变点检测 | Codeboy's Blog

一台机器原来每件产品有 $1/3$ 的概率出现某种标记。机器出问题后，这个概率会升到 $2/3$ 。我们按顺序检查产品，希望尽快发现改变，却不知道它会从哪一件开始。如果把从第一件到现在的全部证据一直相乘，前面很长一段正常记录可能压住刚刚出现的异常。CUSUM 的办法是：同时考虑每一个可能的改变起点，只保留其中最有力的一段证据。它看起来需要保存许多候选段，实际上只需更新一个数。

这个方法还带来一个容易误读的数：正常运行时，平均等多久才误报。平均等待很长，并不表示永远误报的概率很小。本页从同一个三状态例子算出这两种量，说明阈值究竟保证什么。

形式陈述

未知的是起点，改变前后的模型已经给定

令 $X_1,X_2,\ldots$ 相互独立， $\mathcal F_n=\sigma(X_1,\ldots,X_n)$ 。改变前后分布分别有共同支配密度 $f_0,f_1$ 。 $\mathbb P_\infty$ 表示始终使用 $f_0$ ； $\mathbb P_\nu$ 表示 $X_1,\ldots,X_{\nu-1}$ 使用 $f_0$ ，从第 $\nu$ 个观测开始使用 $f_1$ 。所以本页的 $\nu$ 是第一个改变后的下标。

先假定两种分布相互绝对连续，且 $f_0\ne f_1$ 。单步似然比和对数增量为

L_n=\frac{f_1(X_n)}{f_0(X_n)},\qquad Z_n=\log L_n. \tag{1}

在时刻 $n$ ，假设改变发生在 $k\le n$ ，相对于始终正常模型的似然比是 $L_{k:n}=\prod_{i=k}^n L_i$ 。取所有候选起点的最大值：

V_n=\max_{1\le k\le n}L_{k:n},\qquad V_0=1,\qquad V_n=\max(1,V_{n-1})L_n. \tag{2}

递推的证明只需分两类：起点 $k=n$ 给出 $L_n$ ，较早起点的最佳值是 $V_{n-1}L_n$ 。新的候选段总可以从当前观测重新开始。

常用的反射对数形式是

C_0=0,\qquad C_n=\max(0,C_{n-1}+Z_n). \tag{3}

归纳可得 $e^{C_n}=\max(1,V_n)$ 。因此，对 $A>1$ ，

\tau_A=\inf\{n\ge1:C_n\ge\log A\} =\inf\{n\ge1:V_n\ge A\} \tag{4}

是同一个停时。这里取“达到或超过”，并约定空集的下确界为 $\infty$ 。 $V_n$ 可以小于一，而 $e^{C_n}$ 总至少为一；不要把两者的数值直接写成相等。若允许 $L_n=0$ ，式 (2) 仍有意义；式 (3) 要以 $\log0=-\infty$ 理解，而不应让程序计算普通浮点对数零。

反射是把过去最低的累计和扣掉

设 $S_0=0$ 、 $S_n=\sum_{i=1}^nZ_i$ 。所有后缀和满足

C_n=S_n-\min_{0\le j\le n}S_j. \tag{5}

证明是将 $\sum_{i=k}^n Z_i=S_n-S_{k-1}$ 取最大值，并加入空段的零。这个公式也解释了重置：当新的 $S_n$ 成为历史最低点时， $C_n=0$ ，此前的不利证据不再拖累未来候选段。

在报警时若令 $j_*$ 为 $S_0,\ldots,S_{\tau_A-1}$ 的最后一个最小点， $j_*+1$ 是一个最大似然的候选改变起点。选择最后一个只是解决并列的约定。它不是对真实改变时刻的置信区间，也没有自动给出定位误差概率。

用有限状态转移，精确算出等待时间

若递推在报警前只有 $m$ 个可能状态，可把它看作Markov 链。令 $Q_{ij}$ 为一步从暂态 $i$ 到暂态 $j$ 的概率；一行缺少的概率就是这一步报警的概率。设 $\tau$ 是从当前状态起，还要等待的观测数， $\mathbf1$ 为全一列向量，则

\mathbb P_i(\tau>n)=e_i^\top Q^n\mathbf1,\qquad u_i=\mathbb E_i\tau. \tag{6}

第一式来自连续 $n$ 步都留在暂态。这里 $n=0$ 时右端是一，符合下一次检查至少还需一步。

如果存在整数 $b$ 和 $\varepsilon>0$ ，使每个暂态在接下来的 $b$ 步内都有至少 $\varepsilon$ 的报警概率，那么

\|Q^{kb}\|_\infty\le(1-\varepsilon)^k. \tag{7}

把时间按块分开，每块仍存活的条件概率至多 $1-\varepsilon$ ，即得这个界。因此 $\sum_{n\ge0}Q^n$ 收敛，望远镜相消给出它等于 $(I-Q)^{-1}$ 。对整数非负等待时间作尾概率求和，得到

\boldsymbol u=(I-Q)^{-1}\mathbf1, \qquad \boldsymbol u=\mathbf1+Q\boldsymbol u. \tag{8}

所以“解一个线性方程”有明确的前提。若根本无法从某些状态报警，矩阵可能不可逆；不能把一个形式逆矩阵当作有限均值证明。

直觉

固定起点的SPRT比较的是“从开始到现在，哪一种固定模型更像”。变点检测问的是“最近是否有一段已经变了”。CUSUM 的最大后缀允许把前面正常的观测放在候选起点之前，但它仍然使用事先指定的 $f_0,f_1$ 。若改变后的参数也未知，还需要另一个估计或混合步骤，本页没有把它藏进公式中。

本页的 $\log$ 与下面两项 KL 散度统一使用自然对数；链接页默认的底二版本须乘以 $\log2$ 才得到这里的单位。在对数可积的条件下，KL 散度给出

\mathbb E_\infty Z_1=-D(f_0\|f_1)<0, \qquad \mathbb E_1 Z_1=D(f_1\|f_0)>0. \tag{9}

正常时累计证据通常往下走，反射把它留在零附近；改变后通常往上走。不过，“通常向下”不等于“永远到不了上边界”。反射不断提供新机会，偶尔连续出现一串有利观测，仍会触发报警。

图左的折线展示一条观测路径，不是均值轨迹。图右把相同递推压缩成三个暂态，因而能算所有路径的概率。图中的回到零也只丢掉当前候选段的累计证据，没有声称机器被修好了。

例子与边界

三个状态，两个不同的运行长度

回到开头的 Bernoulli 例子： $f_0$ 的成功概率为 $1/3$ ， $f_1$ 为 $2/3$ 。成功时 $L=2$ ，失败时 $L=1/2$ 。用 $\log2$ 作单位，状态每次上升一格或下降一格，下降到零后不再变负。

取 $A=8$ ，即第三格报警。记当前实际成功概率为 $p$ 、 $q=1-p$ ，暂态按 $0,1,2$ 排列：

Q(p)=\begin{pmatrix}q&p&0\\q&0&p\\0&q&0\end{pmatrix}. \tag{10}

例如从第二格失败，回到第一格；从第二格成功，直接报警，所以第三行之和为 $q$ 。均值方程逐行是

u_0=1+qu_0+pu_1,\quad u_1=1+qu_0+pu_2,\quad u_2=1+qu_1. \tag{11}

任意暂态接下来连续三次成功都会报警，概率至少 $p^3>0$ ，因而前面的有限性证明适用。精确求解得到

\boldsymbol u(1/3)=(33,30,21)^\top, \qquad \boldsymbol u(2/3)=\tfrac18(51,39,21)^\top. \tag{12}

从零开始，正常时平均 $33$ 次检查才报警；若一开始就已经改变，平均检测延迟是 $51/8=6.375$ 次。改变发生前若累计状态已在第一格，之后的平均延迟是 $39/8$ ，不是重新从零计算的 $51/8$ 。

最坏历史延迟怎样得到

在给定同一串改变后观测的情况下，从较高初值出发的反射过程始终不低于从零出发的过程。这是映射 $c\mapsto\max(0,c+z)$ 的单调性逐步传递的结果。因此较高初值的报警不晚于零初值。

对这个例子，在 $\tau_A\ge\nu$ 的历史上，改变前状态是 $0,1,2$ 之一；独立性使改变后的未来不再依赖更早历史。已经提前报警的历史给出零延迟。于是

\sup_{\nu\ge1}\operatorname*{ess\,sup} \mathbb E_\nu[(\tau_A-\nu+1)^+\mid\mathcal F_{\nu-1}] =\frac{51}{8}. \tag{13}

上界来自单调耦合；等号在 $\nu=1$ 、确定的零初值处取得。这计算了当前规则的最坏历史延迟，并没有证明它优于所有其他规则。CUSUM 的更强最优性定理需要另述准则和条件，本页不借用那个结论来替代校准。

平均误报等待，不是永远不误报的概率

正常时连续三次成功的概率是 $1/27$ 。不重叠的三次观测块相互独立；前 $k$ 块都不是全成功的概率为 $(26/27)^k$ 。只要有一块全成功，不论块前状态如何，都已在块内报警。所以

\mathbb P_\infty(\tau_8<\infty)=1, \qquad \mathbb E_\infty\tau_8=33. \tag{14}

两式完全相容。不能把阈值 $8$ 解读为“无限期误报概率至多 $1/8$ ”。如果实际任务只监测十次，应算相应的有限时间量：

\mathbb P_\infty(\tau_8\le10) =1-e_0^\top Q(1/3)^{10}\mathbf1 =\frac{13609}{59049}. \tag{15}

同样的现象不限于三状态链。若两个独立同分布模型不同，且 $\mathbb E_\infty L=1$ ，则存在 $c>1$ 满足 $\mathbb P_\infty(L\ge c)>0$ ；否则 $L\le1$ 且均值一会迫使 $L=1$ 几乎处处。足够长的全有利块使后缀积达到任意固定 $A$ ，并给出几何块尾界。此事件在改变后也有正概率，因为其概率是 $\mathbb E_\infty[L\mathbf1_{\{L\ge c\}}]$ 。因此两种模型下都能得到有限均值。相反，若模型完全相同， $L\equiv1$ ，CUSUM 在 $A>1$ 时永不报警。

推论与应用

与 SR 的比较给出一个保守预算

Shiryaev–Roberts 检测把相同的候选后缀积求和，而非取最大。零初值时记其统计量为 $R_n$ ，则逐路径有 $R_n\ge V_n$ 。对相同 $A>1$ ，它的报警时刻 $\tau_A^{\rm SR}\le\tau_A$ 。由 SR 页的有界停止论证，正确正常模型下

\mathbb E_\infty\tau_A\ge \mathbb E_\infty\tau_A^{\rm SR}\ge A, \qquad \mathbb P_\infty(\tau_A\le N)\le\min(1,N/A). \tag{16}

后一个界来自 $A\mathbb P(\tau_A^{\rm SR}\le N)\le\mathbb E(\tau_A^{\rm SR}\wedge N)\le N$ 。它不是一个对无限期有用的固定小概率界。对三状态例子，精确的 $33$ 比保守下界 $8$ 更有校准价值。

最大和求和也会产生真正不同的规则。在公平硬币正常模型、改变后必定成功的模型中，成功给 $L=2$ ，失败给零。阈值 $A=6$ 下，CUSUM 要连续三次成功，平均等待 $14$ 次；SR 两次成功已有 $R=6$ ，平均只等 $6$ 次。这个比较说明相同数值阈值不是相同的误报约束，不能据此宣布哪个方法在公平条件下更快。

实施时需要交付什么

在线递推每个观测只用常数次算术操作与一个当前状态；若还需候选起点，只需额外记住最近重置位置。有限状态校准则是另一个计算： $m$ 状态的稠密线性求解约需 $O(m^3)$ 次域运算，递推 $N$ 步存活分布约需 $O(Nm^2)$ 次；本例矩阵稀疏，可以更快。精确有理数的分子分母还会增长，不能把这些操作数当作位复杂度。

运行长度证书终点要求同时交付：递推与显式后缀的一致性、转移矩阵、有限吸收理由、精确均值、指定时域的报警概率，以及改变初值或阈值后的重算。程序不把一个模拟平均数当作误报保证，也不对未知或失配的 $f_0$ 自动担保。

参考资料

[1] G. V. Moustakides, A. S. Polunchenko, A. G. Tartakovsky, Numerical Comparison of CUSUM and Shiryaev–Roberts Procedures for Detecting Changes in Distributions, 2009 作者版本，§§2.1–2.2 与 §3 的首次转移方程。本文独立推导离散有限状态例子，不采用其连续模型的数值近似误差结论。

[2] G. V. Moustakides, Optimal Stopping Times for Detecting Changes in Distributions, Annals of Statistics 14 (1986), 1379–1387。核读 pp.1379–1380 的模型、最坏历史延迟和单调性。本文只计算给定规则的延迟，不援引该文后面的最优性证明，也不把其连续似然条件用于离散例子。

形式陈述 ​

未知的是起点，改变前后的模型已经给定 ​

反射是把过去最低的累计和扣掉 ​

用有限状态转移，精确算出等待时间 ​

三个状态，两个不同的运行长度 ​

最坏历史延迟怎样得到 ​

平均误报等待，不是永远不误报的概率 ​

与 SR 的比较给出一个保守预算 ​

实施时需要交付什么 ​

反向链接