当前位置：首页 > web >正文

方差计算算法-在线更新算法

web 2025/7/3 20:17:18

from: https://en.m.wikipedia.org/wiki/Algorithms_for_calculating_variance#

方差计算算法在计算统计学中起着重要作用。设计一个好的方差计算算法的关键困难是，方差公式可能涉及平方和，这可能导致数值不稳定，以及在处理大值时出现算术溢出。（p.s. 数值稳定性）

常规算法

计算大小为N的数据的方差：

$eq?%5Csigma%20%5E%7B2%7D%20%3D%20%5Cfrac%7B%5Csum%7B%28X%20-%20%5Cmu%20%29%5E2%7D%7D%7BN%7D$ , (1.1)

其中 $eq?%5Csigma%20%5E%7B2%7D$ 为总体方差, $eq?X$ 为数据， $eq?%5Cmu$ 为均值， $eq?N$ 为数据多少。

由 $eq?D%28X%29%20%3D%20E%28%28X%20-%20E%28X%29%29%5E2%29%20%3D%20E%28X%5E2%29%20-%20E%5E2%28X%29$ 可以推得

$eq?%5Csigma%20%5E2%20%3D%20%5Cbar%7Bx%5E2%7D%20-%20%28%5Cbar%7Bx%7D%29%5E2%20%3D%5Cfrac%7B%5Csum_%7Bi%3D1%7D%5EN%20x_%7Bi%7D%5E%7B2%7D%20-%20%5Cfrac%7B%28%5Csum_%7Bi%3D1%7D%5EN%20x_%7Bi%7D%29%5E2%7D%7BN%7D%7D%7BN%7D$ (1.2)

使用贝塞尔校正来计算n个数据的有限样本的总体方差的无偏估计，公式为：

$eq?s%20%5E2%20%3D%28%5Cfrac%7B%5Csum_%7Bi%3D1%7D%5En%20x_%7Bi%7D%5E%7B2%7D%7D%7Bn%7D%20-%28%5Cfrac%7B%5Csum_%7Bi%3D1%7D%5En%20x_%7Bi%7D%7D%7Bn%7D%29%5E2%20%29%5Ccdot%20%5Cfrac%7Bn%7D%7Bn-1%7D$ (1.3)

(相当于1.2中的分母的N替换为n-1,分子的N替换为n)

总结：

$%28n-1%29$

此算法可以适用于有限样本数的总体方差，只需要把分母的n-1替换成n。但是由于 $eq?SumSq$ 和 $n$ 可能会非常接近，计算过程中的取舍可能会导致精度问题，因此不推荐在实践过程中使用此算法。

计算偏移数据

方差对于位置参数的变化不会产生变化，利用此特性可以避免上述计算过程中的舍入问题：

根据 $eq?Var%28X-K%29%3DVar%28X%29$ p.s. 其中的K可以是任何常数

可以得到新的方差计算公式：

$eq?%5Csigma%20%5E2%20%3D%20%5Cfrac%7B%5Csum_%7Bi%3D1%7D%5En%28x_%7Bi%7D-K%29%5E%7B2%7D%20-%20%5Cfrac%7B%28%5Csum_%7Bi%3D1%7D%5En%28x_%7Bi%7D-K%29%29%5E2%7D%7Bn%7D%7D%7Bn-1%7D$ (1.4)

(1.4)中的K越接近均值，则方差的计算结果越精确；实际计算过程中使用数据中的任意值都能够保证方差计算的数值稳定性。

如果 $eq?%28x_%7Bi%7D-K%29$ 的值很小，则其平方和没有问题;相反，如果值很大很大，则必然意味着方差也很大。在任何情况下，公式中的第二项总是小于第一项，因此不可能发生抵消。

算法中的K值设定为数据中的第一项，下面给出python代码：

def shifted_data_variance(data):if len(data) < 2:return 0.0K = data[0]n = Ex = Ex2 = 0.0for x in data:n += 1Ex += x - KEx2 += (x - K) ** 2variance = (Ex2 - Ex**2 / n) / (n - 1)# use n instead of (n-1) if want to compute the exact variance of the given data# use (n-1) if data are samples of a larger populationreturn variance

p.s. 计算数据精确值时把 $eq?%28n-1%29$ 替换为 $eq?n$ ;计算大数据量样本时使用 $eq?%28n-1%29$ 。

此方法也可用于增量计算（指一部分数据发生变化时，只计算变化部分的数据而无需重复计算全部数据），下面给出python代码：

##############################################################
## 这里是参数设置
##执行增量计算需要设置好位置参数K、均值Ex、Ex2和已经计算的数据个数
##默认情况下参数为0，此时表示还未计算方差
K = Ex = Ex2 = 0.0   
n = 0
##############################################################def add_variable(x):   #增加数据（单个）global K, n, Ex, Ex2if n == 0:K = xn += 1Ex += x - KEx2 += (x - K) ** 2def remove_variable(x):    #减少数据（单个）global K, n, Ex, Ex2n -= 1Ex -= x - KEx2 -= (x - K) ** 2def get_mean():global K, n, Exreturn K + Ex / ndef get_variance():global n, Ex, Ex2return (Ex2 - Ex**2 / n) / (n - 1)
################################################################执行方式：    
## example.    data = [1,2,3,4,5,6]
##1.依次加入data中的数据
data = [1,2,3,4,5,6]
for x in data: add_variable(x)print( get_mean() )
print( get_variance() )##2.假设已经使用shifted_data_variance计算出了data中前面5个元素的方差
##  K = 1 ； Ex = 10.0 ；Ex2 = 30.0 ；  n = 5
K,Ex,Ex2,n = (1,10.0,30.0,5)
add_variable(6)   #增量计算，加入第6个元素print( get_mean() )
print( get_variance() )

两步算法

一个替代方案是，使用不同的方差公式。同样的，第一步计算样本均值：

$\bar{x}=\frac{\sum {^n_{j=1}x_{j}}}{n}$

第二步，根据均值计算差的平方和：

$sample\,variance=s^{2}=\frac{\sum {_{i=1}^n}(x_{i}-\bar{x})^{2}}{n-1}$ , 其中s表示标准差。python代码：

def two_pass_variance(data):n = len(data)mean = sum(data) / nvariance = sum([(x - mean) ** 2 for x in data]) / (n - 1)return variance

Welford 在线算法

此算法计算时，样本数据 $x_i$ 只需要使用一次，这样可以减少数据存储。算法的主要思想是建立起递推关系式。(下面均为计算第n个样本数据时的均值和方差)

均值： $\bar{x_n}=\bar{x_{n-1}}+\frac{x_n-\bar{x_{n-1}}}{n}$

有偏方差： $\sigma _{n}^{2} = \sigma _{n-1}^{2} + \frac{(x_n - \bar{x_{n-1}})(x_n - \bar{x_{n}}) - \sigma _{n-1}^{2}}{n}$

无偏估计： $s_{n}^{2} = s_{n-1}^{2} + \frac{(x_n-\bar{x_{n-1}} )^2}{n} -\frac{s_{n-1}^{2}}{n-1}$

需要注意的是上面的计算方法存在数值不稳定问题，记 $M_{2,n} = \sum_{i=1}^{n} (x_i-\bar{x_n})^2$ 更好的计算方法如下：

均值： $M_{2,n} = M_{2,n-1}+(x_n-\bar{x_{n-1}} )(x_n-\bar{x_{n}})$

有偏方差： $\delta _n^2 = \frac{M_{2,n}}{n}$

无偏估计： $s_n^2 = \frac{M_{2,n}}{n-1}$

公式推导：

（均值推导略）由方差 $D(X)=E((X-E(X))^2)=E(X^2)-E^2(X)$ 可知，

$\sigma_n ^{2}=\frac{1}{n} \sum_{i=1}^{n} x_i^2-(\bar{x_n} )^2$ 和 $\sigma_{n+1} ^{2}=\frac{1}{n+1} \sum_{i=1}^{n+1} x_i^2-(\bar{x_{n+1}} )^2$ 。从而得到

$n\sigma_n ^{2}=\sum_{i=1}^{n} x_i^2-n(\bar{x_n} )^2$ 和 $(n+1)\sigma_{n+1} ^{2}= \sum_{i=1}^{n+1} x_i^2-(n+1)(\bar{x_{n+1}} )^2$ ，两个结合可以得到 $(n+1)\sigma_{n+1} ^{2}= n\sigma_n ^{2}+n(\bar{x_n} )^2 + x_{n+1}^2-(n+1)(\bar{x_{n+1}} )^2$ ,