迭代递推计算均值、方差的无偏估计（含C++实现）

前言

对于一个序列而言，求均值和方差根据定义式是不难的，其时空复杂度均为 $\mathcal{O}\left(N\right)$ 。但有的时候，我们的样本是一个一个给的，此时新来了一个样本，我们总不可能把原来的样本都捞出来再算一次均值、方差吧，那样时空复杂度都是 $\mathcal{O}\left(N\right)$ 了。因此，我们需要一个递推的方式，假设我们已知前 $n$ 个样本的均值和方差 ${\hat{\mu }}_n,{\hat{\sigma }}_n^2$ ，且知道了新的样本 $x_{n+1}$ ，以 $\mathcal{O}\left(1\right)$ 复杂度给出 ${\hat{\mu }}_{n+1},{\hat{\sigma }}_{n+1}^2$ 。

数学推导

我们知道样本均值、样本方差的无偏估计如下

$${\hat{\mu }}_n=\frac{1}{n}\sum _{k=1}^n{x}_k\text{\hspace{1em}(1.1)}$$

$${\hat{\sigma }}_n^2=\frac{1}{n-1}\sum _{k=1}^n{\left(x_k-{\hat{\mu }}_n\right)}^2\text{\hspace{1em}(1.2)}$$

不妨令

$${\beta }_n=\frac{1}{n}\text{\hspace{1em}(1.3)}$$

则均值递推式如下

$$\begin{array}{rl}{\hat{\mu }}_{n+1}& =\frac{1}{n+1}\sum _{k=1}^{n+1}{x}_k\\ & =\frac{n}{n+1}\frac{1}{n}\left(\sum _{k=1}^n{x}_k+x_{n+1}\right)\\ & =\frac{n}{n+1}\frac{1}{n}\sum _{k=1}^n{x}_k+\frac{1}{n+1}{x}_{n+1}\\ & =\left(1-{\beta }_{n+1}\right){\hat{\mu }}_n+{\beta }_{n+1}{x}_{n+1}\\ & ={\hat{\mu }}_n+{\beta }_{n+1}\left(x_{n+1}-{\hat{\mu }}_n\right)\end{array}\text{\hspace{1em}(1.4)}$$

方差递推式的推导要复杂一些，首先将式 $\left(1.4\right)$ 稍微变形，引入式 $\left(1.5\right)$ 。

$$\begin{array}{rl}{x}_{n+1}-{\hat{\mu }}_{n+1}& =\left(x_{n+1}-{\hat{\mu }}_n\right)+\left({\hat{\mu }}_n-{\hat{\mu }}_{n+1}\right)\\ & =\left(x_{n+1}-{\hat{\mu }}_n\right)-{\beta }_{n+1}\left(x_{n+1}-{\hat{\mu }}_n\right)\\ & =\left(1-{\beta }_{n+1}\right)\left(x_{n+1}-{\hat{\mu }}_n\right)\end{array}\text{\hspace{1em}(1.5)}$$

然后引入式 $\left(1.6\right)$ 。

$$\begin{array}{rl}{\beta }_n{\left(x_{n+1}-{\hat{\mu }}_{n+1}\right)}^2& ={\beta }_n{\left(1-{\beta }_{n+1}\right)}^2{\left(x_{n+1}-{\hat{\mu }}_n\right)}^2\\ & =\frac{1}{n}\frac{n}{n+1}\frac{n}{n+1}{\left(x_{n+1}-{\hat{\mu }}_n\right)}^2\\ & ={\beta }_{n+1}\left(1-{\beta }_{n+1}\right){\left(x_{n+1}-{\hat{\mu }}_n\right)}^2\\ & =\left[{\beta }_{n+1}-{\left({\beta }_{n+1}\right)}^2\right]{\left(x_{n+1}-{\hat{\mu }}_n\right)}^2\end{array}\text{\hspace{1em}(1.6)}$$

然后再对样本方差定义式化简，得到式 $\left(1.7\right)$ 。

$$\begin{array}{rl}{\hat{\sigma }}_{n+1}^2& =\frac{1}{n}\sum _{k=1}^{n+1}{\left(x_k-{\hat{\mu }}_{n+1}\right)}^2\\ & =\frac{1}{n}{\left(x_{n+1}-{\hat{\mu }}_{n+1}\right)}^2+\frac{n-1}{n}\left(\frac{1}{n-1}\sum _{k=1}^n{\left[\left(x_k-{\hat{\mu }}_n\right)+\left({\hat{\mu }}_n-{\hat{\mu }}_{n+1}\right)\right]}^2\right)\\ & ={\beta }_n{\left(x_{n+1}-{\hat{\mu }}_{n+1}\right)}^2+\left(1-{\beta }_n\right)\frac{1}{n-1}\sum _{k=1}^n{\left[\left(x_k-{\hat{\mu }}_n\right)+\left({\hat{\mu }}_n-{\hat{\mu }}_{n+1}\right)\right]}^2\\ & ={\beta }_n{\left(x_{n+1}-{\hat{\mu }}_{n+1}\right)}^2+\left(1-{\beta }_n\right)\frac{1}{n-1}\sum _{k=1}^n{\left(x_k-{\hat{\mu }}_n\right)}^2+{\left({\hat{\mu }}_n-{\hat{\mu }}_{n+1}\right)}^2\\ & =\left[{\beta }_{n+1}-{\left({\beta }_{n+1}\right)}^2\right]{\left(x_{n+1}-{\hat{\mu }}_n\right)}^2+\left(1-{\beta }_n\right){\hat{\sigma }}_n^2+{\left({\beta }_{n+1}\right)}^2{\left(x_{n+1}-{\hat{\mu }}_n\right)}^2\\ & =\left(1-{\beta }_n\right){\hat{\sigma }}_n^2+{\beta }_{n+1}{\left(x_{n+1}-{\hat{\mu }}_n\right)}^2\end{array}\text{\hspace{1em}(1.7)}$$

式 $\left(1.7\right)$ 的推导第三行到第四行的等号是因为完全平方展开的交叉项为 $0$ ，具体推导可见式 $\left(1.8\right)$ 。

$$\begin{array}{rl}\sum _{k=1}^n\left(x_k-{\hat{\mu }}_n\right)\left({\hat{\mu }}_n-{\hat{\mu }}_{n+1}\right)& =\left({\hat{\mu }}_n-{\hat{\mu }}_{n+1}\right)\sum _{k=1}^n\left(x_k-{\hat{\mu }}_n\right)\\ & =\left({\hat{\mu }}_n-{\hat{\mu }}_{n+1}\right)\left(\sum _{k=1}^n{x}_k-\sum _{k=1}^n{x}_k\right)\\ & =0\end{array}\text{\hspace{1em}(1.8)}$$

最终得到递推关系式（化简成这样为了保留 $\left(x_{n+1}-{\hat{\mu }}_n\right)$ 公共项，减少计算量）

$$\{\begin{array}{l}{\hat{\mu }}_{n+1}={\hat{\mu }}_n+{\beta }_{n+1}\left(x_{n+1}-{\hat{\mu }}_n\right)\\ {\hat{\sigma }}_{n+1}^2=\left(1-{\beta }_n\right){\hat{\sigma }}_n^2+{\beta }_{n+1}{\left(x_{n+1}-{\hat{\mu }}_n\right)}^2\end{array}\text{\hspace{1em}(1.9)}$$

C++代码实现

核心代码其实只有这么一点

void SeqMeanVar::AppendImpl(double new_value) {
    double xn1_mun = new_value - m_mean;    // x(n + 1) - mu(n)
    double rev_beta_n = 1 - 1 / m_n;        // 1 - beta(n)
    ++ m_n;
    double beta_n1 = 1 / m_n;               // beta(n + 1)
    m_mean += beta_n1 * xn1_mun;            // mu(n + 1) = mu(n) + beta(n + 1) * (x(n + 1) - mu(n))
    m_var = rev_beta_n * m_var + beta_n1 * xn1_mun * xn1_mun;   // var(n + 1) = (1 - beta(n)) * var(n) + beta(n + 1) * (x(n + 1) - mu(n))^2
}

完整代码（含测试样例）如下——

#include <iostream>     // cout

// 有初始值的均值迭代器
class SeqMeanVar
{
public:
    SeqMeanVar();
    SeqMeanVar(double init_value);
    const double GetN() const;
    const double GetMean() const;
    const double GetVar() const;
    // type=0 为检查均值是否有意义, type=1 为检查方差是否有意义
    bool IsValid(int type=1) const;

    void Append(double new_value);

private:
    void InitialConstruct(double init_value);
    void AppendImpl(double new_value);
    double m_n;
    double m_mean;
    double m_var;

friend std::ostream &operator<<(std::ostream &os, SeqMeanVar &seq);
};

SeqMeanVar::SeqMeanVar()
    : m_n(0)
    , m_mean(0)
    , m_var(0)
{
}

SeqMeanVar::SeqMeanVar(double init_value)
{
    InitialConstruct(init_value);
}

void SeqMeanVar::InitialConstruct(double init_value) {
    m_n = 1;
    m_mean = init_value;
    m_var = 0;
}

const double SeqMeanVar::GetN() const {
    return m_n;
}

const double SeqMeanVar::GetMean() const {
    return m_mean;
}

const double SeqMeanVar::GetVar() const {
    return m_var;
}

bool SeqMeanVar::IsValid(int type) const {
    switch (type) {
        case 0: return m_n >= 1;    // 检查均值
        case 1: return m_n >= 2;    // 检查方差
        default: return false;
    }
}

void SeqMeanVar::Append(double new_value) {
    if (m_n == 0)
        InitialConstruct(new_value);
    else
        AppendImpl(new_value);
}

void SeqMeanVar::AppendImpl(double new_value) {
    double xn1_mun = new_value - m_mean;    // x(n + 1) - mu(n)
    double rev_beta_n = 1 - 1 / m_n;        // 1 - beta(n)
    ++ m_n;
    double beta_n1 = 1 / m_n;               // beta(n + 1)
    m_mean += beta_n1 * xn1_mun;            // mu(n + 1) = mu(n) + beta(n + 1) * (x(n + 1) - mu(n))
    m_var = rev_beta_n * m_var + beta_n1 * xn1_mun * xn1_mun;   // var(n + 1) = (1 - beta(n)) * var(n) + beta(n + 1) * (x(n + 1) - mu(n))^2
}

std::ostream &operator<<(std::ostream &os, SeqMeanVar &seq) {
    os << "(n = " << seq.m_n << ", mean = " << seq.m_mean << ", var = " << seq.m_var << ")";
    return os;
}

int main()
{
    SeqMeanVar seqMV(2);
    std::cout << seqMV << std::endl;
    seqMV.Append (4);
    std::cout << seqMV << std::endl;
    seqMV.Append (6);
    std::cout << seqMV << std::endl;
    seqMV.Append (7);
    std::cout << seqMV << std::endl;
    seqMV.Append (8);
    std::cout << seqMV << std::endl;

    getchar();
    return 0;
}