方差这个东西凡是搞数据的人都绕不开。日常处理一批数据大家习惯用现成工具一把梭均值方差直接出结果。但真碰到数据是流式产生的、或者数据量大到内存装不下、又或者你需要实时更新统计量的时候教科书上那个公式就不太好使了必须换一套思路。这篇文章我想把方差迭代计算公式这件事讲透包括它怎么推导、怎么写成代码、数值稳定性到底有多重要以及我踩过的那些坑。这套方法在我实际工作中救过我好几次尤其是做在线监控指标和增量数据处理那阵子。如果你也写数据管道、做实时特征工程或者在嵌入式设备上做均值方差统计这篇文章应该能帮你少走不少弯路。1. 为什么需要迭代计算方差1.1 教科书公式的痛点两遍遍历和内存占用先复习一下大学统计课上学的东西。给定一组数据 x₁, x₂, ..., xₙ总体方差的无偏估计样本方差是s² Σ(xᵢ - x̄)² / (n - 1)其中 x̄ 是均值。这个公式逻辑清晰是定义也是教科书标配。但如果你真的按这个公式去写工程代码很快就会遇到两个现实问题。第一个问题是两遍遍历。必须先算一遍均值 x̄再遍历一遍求每个点与均值的差的平方和。如果数据已经全部在内存里这倒没什么可要是数据是一个一个到达的——比如传感器读数、股票行情、服务器日志指标——你就必须先把数据缓存起来等凑齐一批再计算。这等于放弃了实时性。第二个问题是内存。数据量大到几百GB甚至TB级别时把所有原始数据都存在内存里就为了算个方差太奢侈了。更合理的方式是边读取数据边计算保持一个很小的常数级状态量而不是缓存全部数据。1.2 迭代公式能解决什么问题迭代计算方差的核心思想是维护一组不断更新的统计量每来一个新数据点就用这个点更新统计量而不是重新计算一遍。整个过程只需要保存当前样本数 n、当前均值 x̄、当前平方差累加和 M2内存占用是 O(1)。这带来几个很实际的好处实时性每来一个点均值方差都立即反映最新状态不需要等待一批数据凑齐。内存可控不管处理多少数据只存三个数值流式处理百万条或亿万条数据没区别。可增量合并多个分区的统计结果可以两两合并这天然适合分布式计算。数值稳定性好这个后面重点展开迭代更新避免了大数吃小数的问题。这套思路还有个名字叫 Welford 算法。它不复杂但其价值被大多数人低估了。2. 方差迭代公式的推导2.1 从均值迭代说起方差迭代推导的第一步是均值的迭代更新。假设已经有 n 个数据的均值 x̄ₙ现在来了第 n1 个数据 x新的均值可以这样算x̄ₙ₊₁ (n · x̄ₙ x) / (n 1)这个公式可以改写成一个更常用的增量形式x̄ₙ₊₁ x̄ₙ (x - x̄ₙ) / (n 1)为什么写成这个形式因为在浮点数计算里大的基准值加一个小的相对偏移比直接做两个大数相加更稳定。后面讲数值稳定性时还会细说。2.2 平方差累加和的增量更新现在到最关键的地方。定义 M2ₙ 为前 n 个数据的平方差累加和M2ₙ Σᵢ₌₁ⁿ (xᵢ - x̄ₙ)²我们要找的是已知 M2ₙ 和 x̄ₙ来了一个新数据 x 后M2ₙ₊₁ 怎么算直接推导如下。设 d₁ x - x̄ₙ也就是新数据点与旧均值的差。新均值 x̄ₙ₊₁ x̄ₙ d₁ / (n1)。那么M2ₙ₊₁ Σᵢ₌₁ⁿ₊₁ (xᵢ - x̄ₙ₊₁)²拆成前 n 个点和新点两部分M2ₙ₊₁ Σᵢ₌₁ⁿ (xᵢ - x̄ₙ₊₁)² (x - x̄ₙ₊₁)²对前 n 个点的部分把每一项做个变换xᵢ - x̄ₙ₊₁ (xᵢ - x̄ₙ) - (x̄ₙ₊₁ - x̄ₙ) (xᵢ - x̄ₙ) - d₁/(n1)展开平方Σᵢ₌₁ⁿ (xᵢ - x̄ₙ₊₁)² Σᵢ₌₁ⁿ [(xᵢ - x̄ₙ)² - 2(xᵢ - x̄ₙ)·d₁/(n1) d₁²/(n1)²]这里用到一个关键性质Σᵢ₌₁ⁿ (xᵢ - x̄ₙ) 0所以交叉项全部消掉。于是Σᵢ₌₁ⁿ (xᵢ - x̄ₙ₊₁)² M2ₙ n·d₁²/(n1)²再看新点那一项(x - x̄ₙ₊₁)² (d₁ - d₁/(n1))² (n·d₁/(n1))² n²·d₁²/(n1)²两部分相加M2ₙ₊₁ M2ₙ n·d₁²/(n1)² n²·d₁²/(n1)² M2ₙ n(n1)·d₁²/(n1)² M2ₙ n·d₁²/(n1)所以M2ₙ₊₁ M2ₙ (x - x̄ₙ)·(x - x̄ₙ₊₁)这个形式特别漂亮。两个括号一个是新点减旧均值一个是新点减新均值乘起来加到 M2 上就行。代码实现时不需要额外开数组一个乘法一个加法就搞定了。2.3 从迭代值还原方差和标准差拿到 n 和 M2 之后方差和标准差就是简单的代数运算总体方差σ² M2 / n样本方差无偏估计s² M2 / (n - 1)标准差s sqrt(M2 / (n - 1))注意分母的 n-1这是贝塞尔校正用于无偏估计。工程上如果数据量很大n 和 n-1 的差别很小无非偏要求时直接用 n 也行。但不管哪种先把 M2 维护好最后一步想怎么除就怎么除。3. 从公式到代码核心实现3.1 最小可用的Python实现把上面的公式变成代码非常直接。我写了一个最精简的版本class OnlineVariance: def __init__(self): self.n 0 self.mean 0.0 self.M2 0.0 def update(self, x): self.n 1 delta x - self.mean self.mean delta / self.n delta2 x - self.mean self.M2 delta * delta2 def variance(self, ddof1): if self.n - ddof 0: return 0.0 return self.M2 / (self.n - ddof) def stddev(self, ddof1): return self.variance(ddof) ** 0.5 def merge(self, other): # 两个统计对象的合并后面单独讲 if other.n 0: return n1, mean1, M2_1 self.n, self.mean, self.M2 n2, mean2, M2_2 other.n, other.mean, other.M2 self.n n1 n2 self.mean mean1 (mean2 - mean1) * n2 / self.n self.M2 M2_1 M2_2 (mean1 - mean2) ** 2 * n1 * n2 / self.n这个代码跑起来什么样给大家一个验证例子import random random.seed(42) ov OnlineVariance() data [random.gauss(100, 15) for _ in range(10000)] for x in data: ov.update(x) import statistics print(手工计算方差:, statistics.variance(data)) print(迭代计算方差:, ov.variance())实测结果两者在小数点后好几位都能对上。这不是魔法就是同一件事的不同表达方式。3.2 带权重版本的迭代公式有些场景下每个数据点的权重不一样。比如时间序列分析里越新的数据权重越高做衰减统计。带权重的 Welford 算法没有那么通用因为权重的定义方式会直接影响推导结果但有一种带指数衰减的实现方式在工程上很常用。我实际用的版本是这样的class ExponentialVariance: def __init__(self, alpha0.1): self.alpha alpha self.mean None self.M2 0.0 self.weight_sum 0.0 def update(self, x): if self.mean is None: self.mean x self.weight_sum 1.0 return self.weight_sum self.alpha * (1 - self.weight_sum) delta x - self.mean self.mean delta * self.alpha # 这里的M2更新用了更接近EWMA的处理方式 diff x - self.mean self.M2 (1 - self.alpha) * (self.M2 self.alpha * delta * diff) def variance(self): return self.M2 / self.weight_sum这个版本借鉴了指数加权移动平均的思量alpha 越大历史数据衰减越快。不过我要提醒一句带权重的方差在工程上有很多不同的口径你要先想清楚业务上到底想要“全量数据的加权方差”还是“指数衰减窗口内的近似方差”。这两个目标不同代码也不同。3.3 滑动窗口版本的实现思路另一个常见需求是滑动窗口方差比如只看最近 N 个数据的方差。滑动窗口的实现有两种思路。第一种是简单的队列缓存法。维护一个长度为 N 的队列每来一个新点先把旧的点移出去再更新统计量。问题在于“移出去”这个操作没有办法通过简单的迭代公式直接做需要额外的数据结构支持。第二种是分段重算法。把窗口分成若干个块每个块单独维护统计量窗口滚动时滑出多少数据就重新计算对应块的统计量。这个在实现上复杂但在数据流很大的时候效率高。我平常用第一种比较多。因为窗口 N 一般不会太大直接用collections.deque和statistics.variance()重算都行。只有当你每个数据点处理耗时很高、且实时性要求很变态的时候才会去优化成块式方案。一个简化实现from collections import deque class SlidingWindowVariance: def __init__(self, window_size): self.window deque(maxlenwindow_size) def update(self, x): self.window.append(x) if len(self.window) 2: return 0.0 # 直接重算窗口内数据 n len(self.window) mean sum(self.window) / n return sum((v - mean) ** 2 for v in self.window) / (n - 1)这个实现每次更新是 O(N)窗口小的时候完全够用。如果窗口达到几万甚至几十万就要考虑分块统计或者双端维护索引那就是另一个话题了。4. 数值稳定性教科书没告诉你的坑4.1 朴素公式为什么会在特定数据上崩掉很多人在用迭代公式之前用的是教科书方差的等价展开式s² Σxᵢ²/n - (Σxᵢ/n)²这个公式数学上跟原始定义等价但浮点数计算中非常危险。问题出在“大数相减”时的灾难性抵消。举个例子假设数据是 [1000000001, 1000000002, 1000000003]。用双精度浮点数存的话Σxᵢ² 大约等于 3 × 10¹⁸平方项约是 10¹⁸ 级别(Σxᵢ/n)² 也是 10¹⁸ 级别两个 10¹⁸ 级别的数相减得到的结果可能是几千或者几万但真实的方差应该是个位数级别。也就是说有效数字被前面的巨大数值吃掉了最终结果完全失真。Welford 算法为什么能避免这个问题因为它在每一步都让数值保持在“与数据本身同数量级”的状态。它维护的 M2 是每个点相对均值的差的平方累加这个累加量一般在数据方差的量级附近而不是数据平方和的量级。在数据均值为 10⁹、波动为 ±1 的场景下Welford 的 M2 在个位数到几十的范围内增长精度丢失非常有限。4.2 实测一次精度对比我用 Python 做一个对比实验import numpy as np data np.array([1e9 0.1 * i for i in range(1000)], dtypenp.float64) # 朴素公式 s1 np.sum(data ** 2) / len(data) - (np.mean(data)) ** 2 # Welford ov OnlineVariance() for x in data: ov.update(x) s2 ov.M2 / ov.n # 对比高精度基准 from decimal import Decimal, getcontext getcontext().prec 50 data_dec [Decimal(str(x)) for x in data] mean_dec sum(data_dec) / len(data_dec) m2_dec sum((x - mean_dec) ** 2 for x in data_dec) s0 m2_dec / len(data_dec) print(高精度基准:, float(s0)) print(朴素公式:, s1) print(Welford算法:, s2)你可以跑一下结果大概率是朴素公式算出一个离谱的负数或者量级完全错误的数字而 Welford 算法给出的结果与高精度基准在有效数字内一致。负数方差这种诡异结果就是朴素公式在数值上失效的铁证。4.3 和Python内置模块对比Python 3.9 之后的statistics.fmean()和statistics.variance()内部其实已经做了均值稳定的处理。但它们的实现拿到单个元素迭代的话标准库并不直接暴露 Welford 的在线接口。人家是“整个序列给你算”你没法让序列边到达边计算。所以我建议在流式场景下自己实现或者封装一个 Welford 类这十几行代码不亏。5. 扩展从方差到协方差、并行计算5.1 协方差的迭代公式跟方差几乎一模一样的思路可以扩展到两个变量的协方差。维护三个状态量n、均值向量两个、M2其实是 Mxy两个变量差值的乘积累加。class OnlineCovariance: def __init__(self): self.n 0 self.mean_x 0.0 self.mean_y 0.0 self.Cxy 0.0 def update(self, x, y): self.n 1 dx x - self.mean_x dy y - self.mean_y self.mean_x dx / self.n self.mean_y dy / self.n self.Cxy dx * (y - self.mean_y) def covariance(self, ddof1): if self.n - ddof 0: return 0.0 return self.Cxy / (self.n - ddof)协方差、相关系数都可以基于这个迭代对象实时算。做在线机器学习特征监控的时候这个类很好用比如监控两个指标之间的相关性有没有发生漂移。5.2 分布式场景统计量的合并另一个极其实用的扩展是并行化。假设你有 100 个数据分片每个分片各自算出一组 (n, mean, M2)怎么合并成全局的统计量答案是用 Chan 等人提出的并行合并公式。两个统计对象 A 和 B合并公式如下n nA nBmean meanA (meanB - meanA) · nB / nM2 M2A M2B (meanA - meanB)² · nA · nB / n这个公式我前面的OnlineVariance.merge()方法里已经写进去了。关键价值在于分布式框架比如 Spark、Flink 或 MapReduce里每个 worker 独立维护统计量最后的 driver 节点只需要两两合并不需要重算一遍全量数据。这个场景在业务报表系统里非常常见。5.3 增量训练场景的使用机器学习里的特征标准化也经常遇到类似问题。训练数据太大不能一次性读入内存想边读边计算每个特征的均值和标准差等统计完成后再做归一化。Welford 算法正是这个场景的标准解法。每个特征单独维护一个OnlineVariance实例数据流式通过时不断更新。全部数据过完之后每个特征的均值方差都在手边了再用于标准化的 transform 即可。我用这个方法处理过上亿条样本的标准化任务效果很稳。相比先把全量数据落盘再算统计量的老办法流程简化了很多而且节约了一次全量读取耗时。6. 常见问题与调试实录6.1 为什么合并后的结果和全量计算不一致很多人觉得合并公式是精确的为什么实测会有细微偏差因为浮点数的运算顺序不同、每个分组内的求和精度不同合并后的统计量与一次性全量计算会有数位级别的误差。这是浮点数本身的特性不是算法错了。判断标准应该是“误差是否在可接受范围”而不是“是否完全相等”。如果你在做跨节点合并结果需要精确一致建议在节点汇总时用高精度运算比如decimal模块或者接受一个容差范围用于断言验证。6.2 样本方差分母选 n 还是 n-1这是个老生常谈的问题但不同场景真的不一样。如果数据代表全量总体比如统计硬盘上所有文件的大小分布用 n 做分母算总体方差即可。如果数据只是样本想要推断总体方差的无偏估计必须用 n-1。深度学习里的 BatchNorm 计算用的是 n因为是对当前 batch 的活跃分布做归一化并不需要无偏估计。流式计算里的监控预警通常也用 n因为状态量是持续累积的不涉及抽样推断。我在工程上默认写variance(ddof1)如果业务方明确说不需要无偏再传 0。宁可配置显式化也不要默认行为让人猜。6.3 处理海量数据时 M2 溢出的担心有人会问数据几亿条M2 会不会溢出一般来说不会因为 M2 的量级受数据方差和数量共同影响而方差如果有限M2 会随着 n 线性增长。float64 能表示的最大值大约是 1.8×10³⁰⁸要想溢出数据方差和规模的乘积得达到这个量级正常业务场景几乎不可能。真正需要担心的反而是数值小的场景比如所有数据都在 0 附近方差极小这时如果用朴素公式很容易被四舍五入误差磨成 0。Welford 算法在这种情况下表现也优于朴素公式。6.4 初始化时的特殊处理当第一个数据点到来时样本数 n0直接套公式会出现除零问题。我之前见过不止一次大家在实现里忘了这个边界条件。稳妥的做法是第一个点只初始化mean xM2 保持 0n 设为 1然后直接返回不参与任何增量计算。这个边界条件看起来非常基础但真实工程里一个只有一两个样本的流式统计器偶尔被调用到不做好保护就会飘出一个 NaN然后整个监控面板全线飘红。细节决定成败。6.5 快速自检方法写完实现之后我习惯用一组很少的数据验证正确性输入 [1, 2, 3, 4, 5]样本方差应该正好是 2.5。输入 [0, 0, 0, 0]方差必须是 0不能出现 NaN 或负数。输入 [1e9, 1e9 1, 1e9 2]方差应该接近 1朴素公式大概率会翻车Welford 应该得到准确值。把数据分成两半分别算两个统计对象再合并和全量计算结果误差应在 1e-10 以内。空序列统计时variance()可以返回 0 或抛出异常但绝不能让后续流程收到 NaN。这五条过了基础实现基本就没问题。7. 实操经验我建议这样选型平心而论迭代计算方差不应该是你手里唯一的统计工具。数据量小、一次性可以全部载入内存的时候直接statistics.variance(data)是最简单的没必要套一个什么迭代算法。但如果你遇到了下面任何一种情况请把 Welford 拿出来数据是流式的不能等全部到达后再计算数据量非常大无法全部放入内存需要在分布式环境下合并统计量需要周期性输出当前累计统计值均值很大、方差很小对数值精度有要求场景判断准确比多写几行代码重要得多。在我做过的那些实时监控项目里Welford 算法不是炫技是真的省事。它维护的只是一个三元组放到 Redis 里保存状态都没问题进程重启了也能从最近的状态继续更新不需要回放历史数据。这一点在长时间运行的服务里非常实用凌晨服务重启早上看监控指标数据不会断档。另外我有个习惯收到新的数据点后先update()隔一段时间把当前统计量落一次快照。万一后面发现计算逻辑写错了还能从最近一个快照回滚重新开始不至于从头再来。这个小习惯帮我节省过不少排查时间。回到本文标题里的四个字迭代计算方差的本质其实是增量维护状态。它不改变方差的定义改变的是数学上等价但计算路径完全不同的更新策略。这个策略一旦掌握你很快就会发现它不只是能算方差协方差的在线估计、卡尔曼滤波里的协方差更新、贝叶斯推理里的后验更新底子都是同一套思想——用当前状态和新观测做一次局部修正而不是全量重算。我个人的体会是数值算法这种“换个计算顺序结果天差地别”的特性比任何经验总结都更能提醒自己代码写出来能跑只是起点算得准才是真功夫。