Welford算法的推导和实现

介绍

该算法是最初由B.P.Welford于1962年提出的计算样本均值和样本方差的算法。
算法如下¹：

初始化 $M_1=x_1，S_1=0$

对于接下来的样本值$x$，使用递推公式

$M_k=M_{k-1}+\frac{x_k-M_{k-1}}{k}$
$S_k=S_{k-1}+(x_k-M_{k-1})(x_k-M_k)$

其中，$2⩽k⩽n$，第$k$个样本方差估计为$s^2=\frac{S_k}{(k-1)}$

推导

样本均值²

通常计算均值的方法是
$$mean=\frac{total}{count}$$
用数学符号表达为：
$$\bar{x}=\frac{1}{n}\sum _{i=1}^n{x}_i$$
但是传统的计算方法有如下两个弊端：

• 当累积的总和很大，在使用浮点型数据类型时可能造成精度缺失和溢出问题
• 计算时必须保存所有数据

这两个问题可以使用一种增加式方法来解决，当有新值出现，再调整均值和方差。

均值可以表示为前n-1个值的均值加上最后的新值除以n：
$$\bar{x}=\frac{1}{n}(\sum _{i=1}^{n-1}{x}_i+x_n)$$
其中前n-1个值的均值可以表示为total/sum的形式：
$${\bar{x}}_{n-1}=\frac{{\sum }_{i=1}^{n-1}{x}_i}{n-1}$$
变换该公式：
$$\sum _{i=1}^{n-1}{x}_i=({\bar{x}}_{n-1})(n-1)$$
使用上式表示前n-1个值的均值，带入均值计算公式：
$${\bar{x}}_n=\frac{1}{n}((n-1){\bar{x}}_{n-1}+x_n)$$
展开上式：
$${\bar{x}}_n=\frac{n{\bar{x}}_{n-1}-{\bar{x}}_{n-1}+x_n}{n}$$
变换一下：
$${\bar{x}}_n={\bar{x}}_{n-1}+\frac{x_n-{\bar{x}}_{n-1}}{n}$$
所以当新增一个值时，新的均值等于旧的均值再加上$\frac{x_n-{\bar{x}}_{n-1}}{n}$

在得到结论之前，我们可以再推导出一个将要用于方差计算中的恒等式。

我们当前的结论变换一下：
$${\bar{x}}_n-{\bar{x}}_{n-1}=\frac{x_n-{\bar{x}}_{n-1}}{n}$$
两端同时乘以n：
$$n({\bar{x}}_n-{\bar{x}}_{n-1})=x_n-{\bar{x}}_{n-1}$$
再同时乘以-1：
$$n({\bar{x}}_{n-1}-{\bar{x}}_n)={\bar{x}}_{n-1}-x_n$$

样本方差计算³

这节将会介绍如何递增地计算样本方差和标准差。
同样地，这种方法可以不需要保存所有数据，并且是“数值稳定”的，即不会有精度和溢出问题。

样本方差算法的推导比样本均值多一些步骤，但同样容易理解。

通常来说，样本方差的数学公式如下：
$$s^2=\frac{{\sum }_{i=1}^n(x_i-{\bar{x}}_n{)}^2}{(n-1)}$$
两遍同乘以$(n-1)$：
$$(n-1)s^2=\sum _{i=1}^n(x_i-{\bar{x}}_n{)}^2$$
令$d_n^2=(n-1)s^2$：
$$d_n^2=\sum _{i=1}^n(x_i-{\bar{x}}_n{)}^2$$
展开等式右边：
$$d_n^2=\sum _{i=1}^n(x_i^2-2x_i{\bar{x}}_n+{\bar{x}}_n^2)$$
进一步展开等式右边：
$$d_n^2=\sum _{i=1}^n{x}_i^2-\sum _{i=1}^{n}2x_i{\bar{x}}_n+\sum _{i=1}^n{\bar{x}}_n^2$$
提出常数：
$$d_n^2=\sum _{i=1}^n{x}_i^2-2{\bar{x}}_n\sum _{i=1}^n{x}_i+{\bar{x}}_n^2\sum _{i=1}^{n}1$$
又因为：
$$\sum _{i=1}^n{x}_i=n{\bar{x}}_n$$
代入得：
$$d_n^2=\sum _{i=1}^n{x}_i^2-2n{\bar{x}}_n^2+n{\bar{x}}_n^2$$
整理得到：
$$d_n^2=\sum _{i=1}^n{x}_i^2-n{\bar{x}}_n^2$$
同理可得到：
$$d_{n-1}^2=\sum _{i=1}^{n-1}{x}_i^2-(n-1){\bar{x}}_{n-1}^2$$
上述两式相减得：
$$d_n^2-d_{n-1}^2=\sum _{i=1}^n{x}_i^2-n{\bar{x}}_n^2-(\sum _{i=1}^{n-1}{x}_i^2-(n-1){\bar{x}}_{n-1}^2)$$
整理得：
$$d_n^2-d_{n-1}^2=\sum _{i=1}^n{x}_i^2-n{\bar{x}}_n^2-\sum _{i=1}^{n-1}{x}_i^2+(n-1){\bar{x}}_{n-1}^2$$
同类项相减：
$$d_n^2-d_{n-1}^2=x_n^2-n{\bar{x}}_n^2+(n-1){\bar{x}}_{n-1}^2$$
展开最后两项：
$$d_n^2-d_{n-1}^2=x_n^2-n{\bar{x}}_n^2+n{\bar{x}}_{n-1}^2-{\bar{x}}_{n-1}^2$$
重排顺序：
$$d_n^2-d_{n-1}^2=x_n^2-{\bar{x}}_{n-1}^2-n{\bar{x}}_n^2+n{\bar{x}}_{n-1}^2$$
提取公因式：
$$d_n^2-d_{n-1}^2=x_n^2-{\bar{x}}_{n-1}^2+n({\bar{x}}_{n-1}^2-{\bar{x}}_n^2)$$
由于
$$a^2-b^2=(a-b)(a+b)$$
将其运用在括号中的表达式：
$$d_n^2-d_{n-1}^2=x_n^2-{\bar{x}}_{n-1}^2+n({\bar{x}}_{n-1}-{\bar{x}}_n)({\bar{x}}_{n-1}+{\bar{x}}_n)$$
上一小节我们得到：
$$n({\bar{x}}_{n-1}-{\bar{x}}_n)={\bar{x}}_{n-1}-x_n$$
代入得：
$$d_n^2-d_{n-1}^2=x_n^2-{\bar{x}}_{n-1}^2+({\bar{x}}_{n-1}-x_n)({\bar{x}}_{n-1}+{\bar{x}}_n)$$
展开得：
$$d_n^2-d_{n-1}^2=x_n^2-{\bar{x}}_{n-1}^2+{\bar{x}}_{n-1}^2+{\bar{x}}_{n-1}{\bar{x}}_n-x_n{\bar{x}}_{n-1}-{\bar{x}}_n{x}_n$$
相互抵消后得：
$$d_n^2-d_{n-1}^2=x_n^2+{\bar{x}}_{n-1}{\bar{x}}_n-x_n{\bar{x}}_{n-1}-{\bar{x}}_n{x}_n$$
可知
$$(x-a)(x-b)=x^2-bx-ax+ab$$
所以原式变换为：
$$d_n^2-d_{n-1}^2=(x_n-{\bar{x}}_{n-1})(x_n-{\bar{x}}_n)$$
所以
$$d_n^2=d_{n-1}^2+(x_n-{\bar{x}}_{n-1})(x_n-{\bar{x}}_n)$$

过程有些繁琐，但是结果已经得到啦！

为了得到方差$s^2$，只需将$d_n^2$除以$n$或者$(n-1)$即可。
$$s_n^2=\frac{d_n^2}{n-1}$$

标准差
$$s_n=\sqrt{\frac{d_n^2}{n-1}}$$

程序实现

package chapter2;

import java.util.Scanner;

public class Welford {
    private int num;
    private double mu;
    private double sum;

    /**
     * 添加新值
     * @param newValue
     */
    public void addNewValue(double newValue){
        num++;
        double delta = (newValue-mu);
        mu+=delta/num;
        sum+=delta*(newValue-mu);
    }

    public Welford() {
    }

    /**
     * @return 均值
     */
    public double mean() {
        return mu;
    }

    /**
     * 
     * @return 方差
     */
    public double var() {
        if (num <= 1) return Double.NaN;
        return sum/(num-1);
    }

    /**
     * 
     * @return 标准差
     */
    public double stddev() {
        return Math.sqrt(this.var());
    }

    /**
     * 
     * @return 样本量
     */
    public int count(){
        return num;
    }

    @Override
    public String toString() {
        return
                "count=" + num +
                ", mean=" + mu +
                ", stddev=" + stddev()
                ;
    }

    public static void main(String[] args) {
        Welford w = new Welford();
        Scanner sc = new Scanner(System.in);
        while (sc.hasNext()){
            w.addNewValue(Double.parseDouble(sc.next()));
        }

        System.out.printf("n      = %d\n",w.count());
        System.out.printf("mean   = %.5f\n",w.mean());
        System.out.printf("stddev = %.5f\n",w.stddev());
        System.out.printf("var    = %.5f\n",w.var());
        System.out.print(w);
    }

}

---

1. Accurately computing running variance ↩︎

2. 样本均值推导 ↩︎

3. 样本方差推导 ↩︎