方差算出了负数
方差是「离均差的平方的平均」。它是一堆平方加起来除以个数，在数学上绝无可能是负的。可是每一本统计教材第一章都会给出的那条「方便计算」的公式，能在四个整数上算出 −16384。这一章讲这个坑，讲它的两种修法，以及一个反常识的结论：那个被到处推荐的在线算法，其实不是最准的那个。
四个数：
10 000 000 001 10 000 000 002 10 000 000 003 10 000 000 004
它们的总体方差，心算就能得出：均值是 10¹⁰ + 2.5，离均差是 −1.5, −0.5, +0.5, +1.5，平方和是 2.25+0.25+0.25+2.25 = 5，除以 4 得 1.25。
现在用统计教材上那条「一遍就能算完」的公式：
Var = E[x²] − (E[x])² （平方的均值 减去 均值的平方）
问：用 double 照着这条公式算，得到什么？
第二问：如果换成 10⁸ 而不是 10¹⁰，答案会好一些吗？
把 BIG 扫一遍
同样四个数 BIG+1, BIG+2, BIG+3, BIG+4，真方差恒等于 1.25，只改 BIG：
| BIG | 教科书一遍法E[x²]−E[x]² | 两遍法 | Welford | numpy.var |
|---|---|---|---|---|
| 0 | 1.25 | 1.25 | 1.25 | 1.25 |
| 10⁶ | 1.25 | 1.25 | 1.25 | 1.25 |
| 10⁸ | 2 | 1.25 | 1.25 | 1.25 |
| 10⁹ | 0 | 1.25 | 1.25 | 1.25 |
| 10¹⁰ | −16384 | 1.25 | 1.25 | 1.25 |
BIG = 10⁸ 那一行最阴险：它给出 2。这是个完全合法的方差值，没有任何地方看起来不对——如果这是你的监控指标，你会照单全收。而真值是 1.25，错了 60%。
10⁹ 给出 0（「这批数据没有波动」），10¹⁰ 给出 −16384（这个至少还会让人警觉）。
为什么
BIG = 10¹⁰ 时： Σx² = 4.0000000200000003e+20 ulp ≈ 65536 n·mean² 也是同一个量级 真正要的答案是这两个数的差 ÷ 4 = 1.25 也就是说：⟨要从两个 4×10²⁰ 的数里，减出一个 5 来⟩ 而那一带的刻度是 65536。 5 连一格的万分之一都不到，它在存进 double 的时候就没了。 剩下的全是舍入噪声——噪声的符号是随机的， 所以这条公式给出正数、零、负数都毫不奇怪。
这是第 3 章的灾难性相消，穿了一件统计学的外衣。而这件外衣骗过的人特别多，因为「E[x²]−E[x]²」这条式子在数学上完全正确，而且它确实只需要扫一遍数据。
修法一：两遍法
老老实实按定义算，扫两遍：
// 第一遍：算均值
let s = 0;
for (const x of a) s += x;
const mean = s / a.length;
// 第二遍：算离均差的平方和
let acc = 0;
for (const x of a) { const d = x - mean; acc += d * d; }
const variance = acc / a.length;
关键在 x - mean 这一步。这确实是两个大数相减——但这次是良性相消（第 3 章的术语）：x 和 mean 都是精确或几乎精确的，相减的结果 d 是个小数，而它本身没有被放大过。之后所有运算都在这个小数上进行，全程待在刻度最密的地方。
两遍法的问题只有一个：要扫两遍。数据在磁盘上、在网络流里、或者根本存不下时，这就不成立了。
修法二：Welford 在线法
1962 年 B. P. Welford 给出了一个只扫一遍、又不相消的办法。想法是边走边更新均值，同时累加一个「和当前均值算的」平方和：
let n = 0, mean = 0, M2 = 0;
for (const x of a) {
n += 1;
const d = x - mean; // 用「旧均值」算的偏差
mean += d / n; // 更新均值
M2 += d * (x - mean); // ★ 旧偏差 × 新偏差
}
const variance = M2 / n; // 总体方差；样本方差用 M2 / (n-1)
最后一行是全部的巧妙之处：d 用旧均值算，(x − mean) 用新均值算，两者相乘。可以验证这个乘积恰好等于 M2 需要增加的量（这是个恒等式，一行代数就能推出来）。
为什么它不相消？因为 M2 从头到尾累加的都是小量（偏差的乘积），累加器长到的规模是「方差 × n」，而不是「均值² × n」。这两个规模差了 (mean/σ)² 倍——在上面那个例子里就是 6.4×10¹⁹ 倍。
Welford 还有两个附带的好处，让它在工程上几乎不可替代：
- 流式。只要
(n, mean, M2)三个数，可以随时报告当前方差，不需要保存数据。 - 可合并。两个分片各自算完，可以用 Chan–Golub–LeVeque 的合并公式拼起来：
M2 = M2ₐ + M2ᵦ + δ²·nₐ·nᵦ/n。这就是分布式系统里算方差的标准做法。
一个反常识：两遍法其实比 Welford 准
Welford 被推荐得太多，以至于常被当成「最准的那个」。实测一下——十万个来自 N(10⁹, 1) 的样本，真方差（用精确求和算的）是 0.9907025774583511：
| 算法 | 算出来的方差 | 相对误差 |
|---|---|---|
| ① 教科书一遍法 | 9088 | 错了 9000 倍 |
| ② 两遍法 | 0.9907025774694592 | 1.1 × 10⁻¹¹ |
| ③ Welford 在线法 | 0.9907025933335 | 1.6 × 10⁻⁸ |
| ④ 精确求和 | 0.9907025774583511 | — |
两遍法比 Welford 准了三个数量级。原因也不难理解：两遍法用的是最终均值，每一个偏差都是「相对真均值」的；Welford 用的是路过时的均值，前几个样本的偏差是相对一个很不准的临时均值算的，那点误差会一直带到最后。
三条算法各有各的位置，选哪条取决于你能扫几遍：
- 数据能放进内存、能扫两遍 ⇒ 两遍法。最准，最简单。
- 只能扫一遍（流、超大数据、在线统计） ⇒ Welford。够准，可合并。
- 教科书那条
E[x²]−E[x]²⇒ 一个位置都没有。它唯一的优点（一遍）Welford 也有，而且不出负数。
顺带一个免费的补丁：如果因为某些原因必须用一遍法，先减去一个「大概的均值」（第一个样本就行）。Var(x) = Var(x − K)，方差平移不变。这一行就能把 10¹⁰ 那个数量级消掉，问题基本解决。这招叫平移法（shifted data algorithm）。
数值不稳定在统计软件的语境里几乎总是指这一章这件事。但要分清两种不同的「不稳定」：
- 算法不稳定（这一章）：数据本身没问题，是公式选错了。换公式就好，免费。
- 估计量不稳定（统计意义）：样本量太小、方差本身很大、存在离群点。换公式没用，要换数据或换估计量。
看到「方差算出来是负的」时，你可以确定是前者——因为后者不可能产生负数。负方差是一个非常干净的信号：它百分之百意味着你的公式里有相消。
负方差是这本书里最容易复现的现象，四行就够：
import statistics, numpy as np
for BIG in (0, 1e6, 1e8, 1e9, 1e10):
v = [BIG + 1, BIG + 2, BIG + 3, BIG + 4]
s = s2 = 0.0
for x in v: s += x; s2 += x * x
naive = s2 / 4 - (s / 4) ** 2 # 教科书一遍法
mu = s / 4
two = sum((x - mu) ** 2 for x in v) / 4 # 两遍法
print(f'BIG={BIG:>10.0e} 一遍法={naive:>12} 两遍法={two} '
f'numpy={np.var(v)} statistics={statistics.pvariance(v)}')
会打出：
BIG= 0e+00 一遍法= 1.25 两遍法=1.25 numpy=1.25 statistics=1.25 BIG= 1e+06 一遍法= 1.25 两遍法=1.25 numpy=1.25 statistics=1.25 BIG= 1e+08 一遍法= 2.0 两遍法=1.25 numpy=1.25 statistics=1.25 BIG= 1e+09 一遍法= 0.0 两遍法=1.25 numpy=1.25 statistics=1.25 BIG= 1e+10 一遍法= -16384.0 两遍法=1.25 numpy=1.25 statistics=1.25
顺带看一件事：numpy.var 和 Python 的 statistics.pvariance 每一行都对。标准库早就替你把这个坑填了——问题只发生在自己手写公式的时候，而这在写 SQL、写流处理、写监控指标时非常常见。
python3 var.py
没 numpy 也能跑，把那两列删掉即可；statistics 是标准库。在线：python.org/shell。
数据库。老版本的 MySQL、SQLite 以及不少 OLAP 引擎的 VAR() / STDDEV() 用的就是一遍法，在「均值远大于标准差」的列上会给出错的甚至负的结果。PostgreSQL 用的是数值上更稳的做法（内部按 numeric 累加）——这是《当真》那本书里「PG 把类型当真」的又一个例子。遇到可疑的方差，先查一下你的引擎用的是哪条公式。
时间戳、坐标、传感器读数是这个坑的三大产地，因为它们的共同特征就是「一个巨大的基准值 + 一点点变化」。算 GPS 坐标的方差、Unix 时间戳的抖动、压力传感器的噪声——全部踩这个坑。平移法（先减掉第一个样本）在这些场合是一行的救命符。
BatchNorm 的 running variance。深度学习框架里的滑动方差用的都是 Welford 的变体或者带动量的更新，而不是一遍法；如果用 fp16 存这些统计量，第 18 章的问题会叠加进来。这也是为什么 BatchNorm 的统计量在混合精度训练里一律保持 fp32。
金融的波动率。股价方差 = 「围绕一个几百块的均值波动几块钱」，比值不算极端，一遍法通常不会炸。但高频数据（毫秒级价格，均值 40000、波动 0.01）就够呛了。加密货币尤其危险，因为价格量级横跨好几个数量级。
「Var = E[x²] − (E[x])² 是标准公式，教材上写着，用它没问题。」
这条式子在数学上完全正确，在数值上是已知的坏做法。它进教材的原因是历史性的：手算年代，扫两遍数据意味着两次抄表，而 Σx 和 Σx² 可以在一次抄表里同时得到。那是为纸笔优化的公式，不是为浮点数优化的。
更值得警惕的是：它大部分时候是对的。数据均值不大时（BIG ≤ 10⁶）它给出完全正确的答案，测试全过、code review 全过、上线之后跑了两年都没事——直到某天有人把时间戳当特征喂了进来。
判据很简单，一句话：看 mean/σ 的比值。这个比值的平方，就是一遍法会丢掉的位数（以 2 为底）。均值 10¹⁰、标准差 1，比值 10¹⁰，平方 10²⁰ ⇒ 丢掉 66 位 ⇒ double 只有 53 位 ⇒ 一位不剩。
正确答案是 D：−16384。换成 10⁸ 也不对——它给出 2，而且这个错答案更危险，因为它看起来完全正常。
10¹⁰+1 这些数确实被精确存下了。丢失发生在 x*x 那一步：平方之后变成 10²⁰ 量级，ulp 涨到 65536，四个平方数一相加，你要的那个 5 就落进了刻度缝里。「输入是精确的」推不出「中间结果是精确的」——而中间结果决定一切。
B 「差几个 ulp」——差的是 1.25 附近的多少个 ulp？ulp(1.25) ≈ 2.2×10⁻¹⁶，而误差是 16385。换算下来差了 7×10¹⁹ 个 ulp。这个答案的思路是「误差总是很小的」，而这一章正是要说明：相消能让误差和结果完全脱钩，「小误差」这个词失去意义。
C 「0」——这是 BIG = 10⁹ 时的答案，方向对了。猜 C 的人已经意识到「那个 5 会被吃掉」，只是没想到吃剩下的不是零，是噪声。噪声有正有负，所以能出负数。结果是 0 还是 −16384，取决于两个 10²⁰ 的数各自被抹到哪一格——这是纯运气。
这一章的改写有三档，按你的约束选：
var = sum(x*x)/n - (sum(x)/n)**2 —— 一遍，可能给负数
能扫两遍：mu = sum(x)/n; var = sum((x-mu)**2)/n 最准，最好读
只能扫一遍：Welford 三行更新（n, mean, M2），够准，还能流式合并
被现有接口卡住、只能一遍：先减一个常数 K（第一个样本就行），Var(x) = Var(x−K)。加一个减法，问题就没了。
这一招（先减掉一个粗略的基准，再在小量上做计算）在这本书里已经是第三次出现了：第 3 章的 GPS 差分、第 5 章的「存增量不存绝对值」、这一章的平移法。它可能是全书性价比最高的一条改写。
这一章的一句话
方差算出负数不是数据的问题，是公式的问题——一条为纸笔优化的恒等式被搬进了浮点数里；而修法是让计算发生在「偏差」上，而不是发生在「原始值」上。
下一章把同一招搬到 AI 里。softmax([1000, 1001, 1002]) 直接照定义算，得到的是 [NaN, NaN, NaN]——因为 exp(1000) 是 Infinity，而 Infinity / Infinity 是 NaN。这是每一个深度学习框架里都写着、但很多人没注意到的一次改写。而它顺带解释了一个更常见的困惑：为什么 PyTorch 要提供 log_softmax，而不让你写 log(softmax(x))。