一种计算标准差的高效方法：Welford迭代法

Posted 2023-02-14 herr_edoc

标准差（Standard Deviation）是概率统计学上非常重要的概念，它被用来描述一组数据的离散程度。
根据标准差的定义，总体标准差σ表示为：而样本标准差S表示为：
总体标准差还有另外一种表示方式：

以上的两种总体标准差计算方法都有三个明显的问题：

1. 代码实现时需要用一个数组来保存所有要计算的数据，如果数据比较大，那么会耗费非常大的空间。
2. 计算时需要对数组总体进行两次遍历，一次计算平均值，一次计算平方，数据较多时耗时也非常大。
3. 对于数据整体差距不大的情况，由于对接近0的浮点数做了平方计算，在数据非常多时累计损失的精度会较大，最终导致计算错误。

针对以上的3个问题，一个叫Welford的大神提出了一种迭代计算标准差的方法：
首先，将均值和方差初始化为0：M1 = x1，S1 = 0，
然后通过以下迭代公式计算每次来新的数据后的均值和标准差：
Mk = Mk-1+ (xk – Mk-1)/k
Sk = Sk-1 + (xk – Mk-1)*(xk – Mk)。
前k个数据的标准差为：s2 = Sk/(k – 1)。（样本标准差用k-1，总体标准差用k）
该计算方法只需要保存两个数据：Mk 和Sk，耗费的空间非常小。由于使用的是迭代的方式，每次的计算量也是非常小的。
通过MATLAB产生1000 000 000个0-1的浮点数，然后分别采用以上三种方法计算放差，最后可以发现方法三最有，方法一表现也尚可，而方法二损失精度非常多，甚至出现了负数的结果。
以下通过MATLAB计算数据集data（可以自己用随机函数生成）每50个数据的方差，对比Welford法和MATLAB内置的var函数（总体方差，除以N）的计算结果，可以发现两者结果是一致的，大家也可以尝试直接计算一个大数据集的标准差，并对比两种方法的精度。

data = rand(10000, 1);
len = length(data);
Mk = data(1);
Mk_1 = Mk;
Sk = 0;
Sk_1 = Sk;
S = zeros(len, 1);
V = zeros(len, 1);
k = 1;
for i = 1:len
    k = mod(i, 50);
    if k == 0
        k = 50;
    end
    xk = data(i);
    Mk = Mk_1 + (xk - Mk_1)/k;
    Sk = Sk_1 + (xk - Mk_1)*(xk - Mk);
    if mod(i, 50) == 0
        S(i) = Sk / k;
        V(i) = var(data(i-49:i), 1);
        Mk = data(i);
        Sk = 0;
    end
    Mk_1 = Mk;
    Sk_1 = Sk;
end
figure
hold on
grid minor
plot(V)
plot(S)
legend("matlab var", "Welford Method”)

参考资料：百度百科、维基百科