5D Library

第三章 早期预警的局限：信噪比、选择偏差与没有预兆的翻转

5D Library › 复杂系统相变、临界慢化与生态韧性的非线性动力学 · English

第三章 早期预警的局限：信噪比、选择偏差与没有预兆的翻转

两个问题

一个预警方法要有用，必须同时回答两个问题：系统真的接近门槛时，它多常报警（检出率）；系统没有接近门槛时，它多常报警（误报率）。只看第一个问题，任何经常喊「狼来了」的方法都显得很灵。5D Library 第 20 本第六章用自回归序列演示过一条简单的规则（肯德尔 τ>0.5\tau > 0.5 就算预警）；本章在放牧模型上把问题系统化，并讨论两类更根本的困难：不留慢化痕迹的转变，以及只研究已经转变的系统造成的偏差。

ROC 曲线

把肯德尔 τ\tau 当作一个分数，高于某个门槛值就报警。门槛值从高到低移动，检出率与误报率同时上升；把每个门槛值下的（误报率，检出率）画出来，就是 ROC 曲线。曲线下的面积 AUC 等于「随机挑一条真正接近门槛的序列和一条没有接近门槛的序列，前者分数更高」的概率 (Hanley and McNeil 1982)：0.5 相当于抛硬币，1 是完美的区分。一般原理见第 26 本第四章。

证据

用 ROC 评估预警指标。 伯蒂格与黑斯廷斯对一组模拟数据、一组恒化器实验数据与一段冰期记录，分别拟合一个稳定模型与一个接近鞍结分岔的模型，从每个拟合的模型各模拟 500 条序列，用 ROC 曲线比较预警指标：两组的肯德尔 τ\tau 分布严重重叠，一条稳定的模拟序列也能偶然给出 τ=0.7\tau = 0.7。把误报率固定在 5% 时，方差指标对这三组数据的检出率只有 25%、5.0% 与 5.4%，他们基于模型的似然方法为 61%、34% 与 100% (Boettiger and Hastings 2012b)。

推导

放牧模型（σ=0.15\sigma = 0.15，1,000 个时间单位）：「接近」组的压力从 1.0 升到 2.5，停在折叠点 2.604 之前；「对照」组的压力始终为 1.0。每组 200 条序列，用三种方式观测同一过程：每个时间单位一次（1,000 个点）；每 5 个时间单位一次（200 个点）；只记录最后 200 个时间单位（200 个点）。去趋势后计算移动窗口中自相关与方差的肯德尔 τ\tau，报告 AUC 与误报率控制在 5% 时的检出率。

推导

# ROC of Kendall's tau of rolling AC(1), variance (Boettiger & Hastings 2012, J R Soc Interface).
# Grazing model, sigma = 0.15, 1,000 time units; "approaching": c rises 1.0 -> 2.5 (fold at
# 2.604); "null": c stays at 1.0. Three views of the same runs: every unit (1,000 points),
# every 5th unit (200 points), only the last 200 units (200 points). 200 series each.
import math, random
rng = random.Random(6)

def simulate(ramp, T=1000, sigma=0.15, dt=0.1):
    x, out = 8.9, []
    for t in range(T):
        c = 1.0 + (1.5 * t / T if ramp else 0.0)
        for _ in range(10):
            x += ((x * (1 - x / 10) - c * x * x / (x * x + 1)) * dt
                  + sigma * math.sqrt(dt) * rng.gauss(0, 1))
        out.append(x)
    return out

def detrend(y, bw):
    k = [(d, math.exp(-0.5 * (d / bw) ** 2)) for d in range(-2 * bw, 2 * bw + 1)]
    out = []
    for i in range(len(y)):
        ok = [(w, y[i + d]) for d, w in k if 0 <= i + d < len(y)]
        out.append(y[i] - sum(w * u for w, u in ok) / sum(w for w, _ in ok))
    return out

def stats(w):
    m = sum(w) / len(w); v = sum((u - m) ** 2 for u in w)
    return sum((w[i] - m) * (w[i + 1] - m) for i in range(len(w) - 1)) / v, v / len(w)

def kendall(y):
    s = sum((y[j] > y[i]) - (y[j] < y[i]) for i in range(len(y)) for j in range(i + 1, len(y)))
    return s / (len(y) * (len(y) - 1) / 2)

def taus(y):
    r = detrend(y, max(3, len(y) // 40)); win = len(r) // 2; step = max(1, win // 25)
    a, v = zip(*(stats(r[i - win:i]) for i in range(win, len(r) + 1, step)))
    return kendall(a), kendall(v)

def auc(pos, neg):                    # P(score of an approaching series > score of a null one)
    return sum((p > q) + 0.5 * (p == q) for p in pos for q in neg) / (len(pos) * len(neg))

def tpr_at(pos, neg, fpr=0.05):
    thr = sorted(neg)[int(math.ceil((1 - fpr) * len(neg))) - 1]
    return sum(p > thr for p in pos) / len(pos), thr

DESIGNS = (("every unit, 1,000 pts", lambda y: y), ("every 5th unit, 200 pts", lambda y: y[::5]),
           ("last 200 units, 200 pts", lambda y: y[-200:]))

def run(n=200):
    series = {k: [simulate(k) for _ in range(n)] for k in (True, False)}
    return {name: {k: [taus(view(y)) for y in v] for k, v in series.items()}
            for name, view in DESIGNS}

if __name__ == "__main__":
    res = run()
    print("design                    indicator   AUC   detected at 5 % false alarms  (tau cut-off)")
    for name, _ in DESIGNS:
        for j, lab in ((0, "AC(1)"), (1, "variance")):
            pos = [t[j] for t in res[name][True]]; neg = [t[j] for t in res[name][False]]
            tp, thr = tpr_at(pos, neg)
            print(f"{name:25s} {lab:9s}  {auc(pos, neg):.2f}         {tp:4.0%}"
                  f"                    ({thr:+.2f})")
    neg = [t[0] for t in res[DESIGNS[0][0]][False]]
    print(f"null series with AC(1) tau > 0.5: {sum(t > 0.5 for t in neg) / len(neg):.0%}"
          f"; largest null tau {max(neg):+.2f}")
# design                    indicator   AUC   detected at 5 % false alarms  (tau cut-off)
# every unit, 1,000 pts     AC(1)      0.94          66%                    (+0.82)
# every unit, 1,000 pts     variance   0.98          86%                    (+0.79)
# every 5th unit, 200 pts   AC(1)      0.66          16%                    (+0.74)
# every 5th unit, 200 pts   variance   0.90          52%                    (+0.79)
# last 200 units, 200 pts   AC(1)      0.54           5%                    (+0.82)
# last 200 units, 200 pts   variance   0.56           8%                    (+0.77)
# null series with AC(1) tau > 0.5: 26%; largest null tau +0.90
三种观测方式下，用肯德尔 \tau 区分「接近门槛」与「什么都没变」的 ROC 曲线（每组 200 条模拟序列）。对角线相当于抛硬币，竖虚线为 5% 的误报率。

结果有四点。第一，记录密集、跨越整个变化过程时，两个指标都区分得很好：自相关的 AUC 为 0.94，方差为 0.98；在 5% 的误报率下分别检出 66% 与 86%。第二，要把误报率控制在 5%，τ\tau 的门槛要设在约 +0.8，而不是常用的 +0.5：在什么都没变的对照序列中，26% 的自相关 τ\tau 超过 0.5，最大达 +0.90；移动窗口中的统计量高度重叠，τ\tau 的分布比独立数据宽得多。第三，记录的密度与跨度作用不同：同样 200 个点，每 5 个单位记录一次、跨越整个过程时，方差仍然有用（AUC 0.90，检出 52%），自相关大为下降（0.66，检出 16%），因为相邻两点的自相关 e−5λe^{-5\lambda} 本身很小；只记录最后 200 个单位时，两者都接近抛硬币（0.54 与 0.56）。达科斯等人也指出，方差的估计受益于更频繁的观测，自相关的估计需要更长的跨度 (Dakos et al. 2015)。第四，这些数字来自参数已知、噪声恒定、确实存在折叠点的模型，真实数据只会更难。

没有预兆的翻转

证据

噪声诱导的转变。 迪特列夫森与约翰森检验了格陵兰 NGRIP 冰芯中丹斯加德—奥施格尔事件（几十年之内的突然升温）之前是否有方差与自相关的上升：两者都没有显著增加。他们认为这些跳变更可能由噪声诱发，而不是由逐渐接近的分岔引起，期待对它们发出预警可能是「一厢情愿」(Ditlevsen and Johnsen 2010)。黑斯廷斯与怀沙姆用几个生态模型说明，只有有限的一类系统会表现出常用的前导指标，许多非线性生态模型在环境变动下没有光滑的势函数，突变之前不留预兆 (Hastings and Wysham 2010)。

证据

失灵的清单。 达科斯等人列出了预警失灵的多种情形：过程或观测误差太大、驱动变化太快、多个门槛相互作用、监测了错误的变量、非局部分岔；噪声诱导的转变、长瞬态与驱动的大幅跳变根本不经过分岔 (Dakos et al. 2015)。反过来，凯菲等人发现同样的指标也出现在非灾变式的平滑转变之前，所以指标上升只能说明系统对扰动越来越敏感，不能说明即将到来的是哪一种转变 (Kéfi et al. 2013)。

选择偏差：检察官谬误

证据

检察官谬误。 伯蒂格与黑斯廷斯模拟了一个参数恒定、不经过任何分岔、只会偶然崩溃的种群模型，只挑出偶然崩溃了的重复，分析崩溃之前的窗口：其中方差与自相关的肯德尔 τ\tau 高于随机挑选的模拟。他们把这种推理比作检察官谬误，混淆了「有预兆时崩溃的概率」与「崩溃时有预兆的概率」；事后截断数据不容易消除这一偏差，有重复的实验可以避免它 (Boettiger and Hastings 2012a)。

推导

压力 cc 保持不变。A：c=2.3c = 2.3、较强的噪声（σ=0.5\sigma = 0.5），400 条各 2,000 个时间单位的序列，有些会偶然越过分界线而崩溃；比较「崩溃前 400 个时间单位」与「从未崩溃的序列中随机截取的 400 个时间单位」。B：c=2.2c = 2.2、更强的噪声（σ=1.0\sigma = 1.0），观察系统在两个吸引盆之间来回闪烁。

推导

# Transitions without a warning, and warnings without an approach. Grazing model with the
# pressure c held CONSTANT (no bifurcation is approached; fold at 2.604).
# A. Chance (noise-induced) collapses, c = 2.3, sigma = 0.5, 2,000 time units per series.
#    Indicators in the 400 units before each collapse are compared with 400-unit windows
#    from series that did not collapse: selecting on the outcome (Boettiger & Hastings 2012).
# B. Flickering, c = 2.2, sigma = 1.0: the state jumps back and forth between the basins.
import math, random
rng = random.Random(21)

def run(c, sigma, T, x=None, dt=0.1, stop=True):
    x = 6.6 if x is None else x
    out = []
    for t in range(T):
        for _ in range(10):
            x += ((x * (1 - x / 10) - c * x * x / (x * x + 1)) * dt
                  + sigma * math.sqrt(dt) * rng.gauss(0, 1))
            x = max(x, 0.01)
        out.append(x)
        if stop and x < 2.5: break
    return out

def stats(w):
    m = sum(w) / len(w); v = sum((u - m) ** 2 for u in w)
    return sum((w[i] - m) * (w[i + 1] - m) for i in range(len(w) - 1)) / v, v / len(w)

def kendall(y):
    s = sum((y[j] > y[i]) - (y[j] < y[i]) for i in range(len(y)) for j in range(i + 1, len(y)))
    return s / (len(y) * (len(y) - 1) / 2)

def taus(seg, win=200, step=8):          # no detrending needed: c is constant
    a, v = zip(*(stats(seg[i - win:i]) for i in range(win, len(seg) + 1, step)))
    return kendall(a), kendall(v)

if __name__ == "__main__":
    tipped, calm = [], []
    for _ in range(400):
        y = run(2.3, 0.5, 2000)
        if y[-1] < 2.5 and len(y) >= 420:
            tipped.append(taus(y[-402:-2]))     # ends 2 units before the fall
        elif y[-1] >= 2.5:
            s = rng.randrange(0, 1600); calm.append(taus(y[s:s + 400]))
    print(f"A. 400 series: {len(tipped)} collapsed by chance after t >= 420,"
          f" {len(calm)} never collapsed")
    for lab, grp in (("windows before a chance collapse", tipped),
                     ("windows from calm series", calm)):
        ta = [t[0] for t in grp]; tv = [t[1] for t in grp]
        print(f"   {lab:33s} mean tau AC(1) {sum(ta) / len(ta):+.2f},"
              f" variance {sum(tv) / len(tv):+.2f}")
        print(f"   {'':33s} share with tau > 0.5: AC(1) {sum(a > 0.5 for a in ta) / len(ta):.0%},"
              f" variance {sum(b > 0.5 for b in tv) / len(tv):.0%}, both"
              f" {sum(a > 0.5 and b > 0.5 for a, b in grp) / len(grp):.0%}")
    y = run(2.2, 1.0, 4000, stop=False)
    low = [u < 3.0 for u in y]
    flips = sum(1 for u, v in zip(low, low[1:]) if u != v)
    a, v = stats(y)
    a0, v0 = stats(run(2.2, 0.15, 4000, stop=False))
    print(f"B. c = 2.2, sigma = 1.0: {flips} crossings of x = 3 in 4,000 units;"
          f" time in lower basin {sum(low) / len(low):.0%}")
    print(f"   AC(1) {a:.2f}, variance {v:.2f}  (same c, sigma = 0.15: AC(1) {a0:.2f},"
          f" variance {v0:.3f})")
    hist = [0] * 10
    for u in y: hist[min(9, int(u))] += 1
    print("   histogram of x (bins of width 1): " + " ".join(f"{h / len(y):.0%}" for h in hist))
# A. 400 series: 87 collapsed by chance after t >= 420, 286 never collapsed
#    windows before a chance collapse  mean tau AC(1) +0.10, variance +0.18
#                                      share with tau > 0.5: AC(1) 21%, variance 29%, both 17%
#    windows from calm series          mean tau AC(1) +0.01, variance +0.04
#                                      share with tau > 0.5: AC(1) 20%, variance 24%, both 16%
# B. c = 2.2, sigma = 1.0: 140 crossings of x = 3 in 4,000 units; time in lower basin 16%
#    AC(1) 0.91, variance 5.21  (same c, sigma = 0.15: AC(1) 0.68, variance 0.031)
#    histogram of x (bins of width 1): 7% 5% 4% 5% 8% 17% 25% 20% 8% 1%

A：400 条序列中 87 条在第 420 个时间单位之后偶然崩溃，286 条从未崩溃。参数从头到尾没有变，但崩溃前窗口中自相关与方差的 τ\tau 平均为 +0.10 与 +0.18，平静序列中为 +0.01 与 +0.04：只因为挑选了「后来崩溃了」的序列，指标就向上偏了。这个偏差不大（τ\tau 超过 0.5 的比例：自相关 21% 对 20%，方差 29% 对 24%），但方向一致。即使在平静序列中，也有 16% 的窗口两个 τ\tau 都超过 0.5。

B：压力 2.2，离折叠点还远，但在强噪声下，系统在 4,000 个时间单位中 140 次穿越 x=3x = 3，16% 的时间停在下方的吸引盆中；自相关 0.91、方差 5.21，而同样的压力、弱噪声时只有 0.68 与 0.031；直方图有两个峰。指标很「高」，却不是因为接近门槛，而是因为噪声足以把系统反复推过分界线。

证据

闪烁。 达科斯、范内斯与谢弗指出，噪声很强时，系统可能在离分岔还很远时就在两个吸引盆之间闪烁并转变，慢化指标可能失效；闪烁也会使方差、自相关与偏度上升，但它的状态分布可以用来勾画可替代的吸引子 (Dakos et al. 2013)。王荣等人分析了云南洱海沉积物中约 125 年的硅藻记录，报告湖泊转为富营养状态之前约 20 年就出现了闪烁 (Wang et al. 2012)。

真实数据

证据

126 组长期生态数据。 伯思等人用 55 个类群、126 组海洋与淡水长期丰度数据，检验方差与自相关能否预测丰度的非线性转折。真阳性很少：方差 9%（170 例中 16 例），自相关 13%（152 例中 19 例）；假阳性比假阴性更常见（方差 53% 对 38%，自相关 47% 对 40%）；两种指标对同一条序列的判断有 43% 不一致。作者认为这两个指标「一贯不能」预测非线性变化，建议结合多种指标与对人为驱动因素的分析 (Burthe et al. 2016)。达科斯等人 2012 年的方法比较也建议对去趋势带宽与窗口长度做敏感性分析，并结合多种指标 (Dakos et al. 2012)。

注

诚实地使用早期预警。 至少应当：（1）事前写下预警规则与门槛；（2）在大量没有发生转变的序列或时段上运行同样的规则，报告误报率与 ROC；（3）同时看方差与自相关，并检查噪声本身是否在变；（4）对窗口与去趋势带宽做敏感性分析；（5）说明指标不能发现的转变类型。指标上升最好被理解为「系统对扰动越来越敏感，值得查明原因」，而不是「门槛将在某日到来」。第八章的核查表就是这样使用它的。

本章小结

预警指标的价值要用检出率与误报率同时衡量。在放牧模型中，记录密集、跨越整个变化过程时，方差与自相关的肯德尔 τ\tau 区分得很好（AUC 0.98 与 0.94），但要把误报率控制在 5%，τ\tau 的门槛要设在约 0.8；记录变稀时自相关先失效，只看最后一小段时两者都接近抛硬币。噪声诱导的翻转与闪烁可以发生在离门槛很远处；只挑选已经崩溃的序列，会在参数从未改变的系统中「看到」预兆。在 126 组真实数据中，真阳性很少，假阳性比假阴性更常见。下一章换一个角度：不预报门槛何时到来，而问系统离门槛有多远，怎样让这段距离变大。

参考文献

Boettiger, Carl, and Alan Hastings. 2012a. “Early Warning Signals and the Prosecutor’s Fallacy.” Proceedings of the Royal Society B: Biological Sciences 279 (1748): 4734–39. https://doi.org/10.1098/rspb.2012.2085.
Boettiger, Carl, and Alan Hastings. 2012b. “Quantifying Limits to Detection of Early Warning for Critical Transitions.” Journal of the Royal Society Interface 9 (75): 2527–39. https://doi.org/10.1098/rsif.2012.0125.
Burthe, Sarah J., Peter A. Henrys, Eleanor B. Mackay, et al. 2016. “Do Early Warning Indicators Consistently Predict Nonlinear Change in Long-Term Ecological Data?” Journal of Applied Ecology 53 (3): 666–76. https://doi.org/10.1111/1365-2664.12519.
Dakos, Vasilis, Stephen R. Carpenter, William A. Brock, et al. 2012. “Methods for Detecting Early Warnings of Critical Transitions in Time Series Illustrated Using Simulated Ecological Data.” PLoS ONE 7 (7): e41010. https://doi.org/10.1371/journal.pone.0041010.
Dakos, Vasilis, Stephen R. Carpenter, Egbert H. van Nes, and Marten Scheffer. 2015. “Resilience Indicators: Prospects and Limitations for Early Warnings of Regime Shifts.” Philosophical Transactions of the Royal Society B: Biological Sciences 370 (1659): 20130263. https://doi.org/10.1098/rstb.2013.0263.
Dakos, Vasilis, Egbert H. van Nes, and Marten Scheffer. 2013. “Flickering as an Early Warning Signal.” Theoretical Ecology 6 (3): 309–17. https://doi.org/10.1007/s12080-013-0186-4.
Ditlevsen, Peter D., and Sigfus J. Johnsen. 2010. “Tipping Points: Early Warning and Wishful Thinking.” Geophysical Research Letters 37 (19): L19703. https://doi.org/10.1029/2010GL044486.
Hanley, James A., and Barbara J. McNeil. 1982. “The Meaning and Use of the Area Under a Receiver Operating Characteristic (ROC) Curve.” Radiology 143 (1): 29–36. https://doi.org/10.1148/radiology.143.1.7063747.
Hastings, Alan, and Derin B. Wysham. 2010. “Regime Shifts in Ecological Systems Can Occur with No Warning.” Ecology Letters 13 (4): 464–72. https://doi.org/10.1111/j.1461-0248.2010.01439.x.
Kéfi, Sonia, Vasilis Dakos, Marten Scheffer, Egbert H. van Nes, and Max Rietkerk. 2013. “Early Warning Signals Also Precede Non-Catastrophic Transitions.” Oikos 122 (5): 641–48. https://doi.org/10.1111/j.1600-0706.2012.20838.x.
Wang, Rong, John A. Dearing, Peter G. Langdon, et al. 2012. “Flickering Gives Early Warning Signals of a Critical Transition to a Eutrophic Lake State.” Nature 492 (7429): 419–22. https://doi.org/10.1038/nature11655.

Fan, Y. (2026). Phase Transitions, Critical Slowing Down and Resilience in Complex Dynamical Systems: From Saddle-Node Bifurcations and Early-Warning Signals to Cascading Failures, Read Beside Classical Chinese Thought on Incipient Change [复杂系统相变、临界慢化与生态韧性的非线性动力学] (Version 1.0.0). 5D Library. https://doi.org/10.5281/zenodo.23770082 · Licensed under Creative Commons Attribution-NonCommercial-NoDerivatives 4.0 International (CC BY-NC-ND 4.0).