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

• • 前言
  • 数学推导
  • C++代码实现

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

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

μ ^ n = 1 n ∑ k = 1 n x k (1.1) hat{mu}_n=frac{1}{n}sum_{k=1}^n{x_k} tag{1.1} μ^​n​=n1​k=1∑n​xk​(1.1)

σ ^ n 2 = 1 n − 1 ∑ k = 1 n ( x k − μ ^ n ) 2 (1.2) hat{sigma}_{n}^{2}=frac{1}{n-1}sum_{k=1}^{n}{left( x_k-hat{mu}_{n} right) ^2} tag{1.2} σ^n2​=n−11​k=1∑n​(xk​−μ^​n​)2(1.2)

不妨令

β n = 1 n (1.3) beta _n=frac{1}{n} tag{1.3} βn​=n1​(1.3)

则均值递推式如下

μ ^ n + 1 = 1 n + 1 ∑ k = 1 n + 1 x k = n n + 1 1 n ( ∑ k = 1 n x k + x n + 1 ) = n n + 1 1 n ∑ k = 1 n x k + 1 n + 1 x n + 1 = ( 1 − β n + 1 ) μ ^ n + β n + 1 x n + 1 = μ ^ n + β n + 1 ( x n + 1 − μ ^ n ) (1.4) begin{aligned} 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{aligned} tag{1.4} μ^​n+1​​=n+11​k=1∑n+1​xk​=n+1n​n1​(k=1∑n​xk​+xn+1​)=n+1n​n1​k=1∑n​xk​+n+11​xn+1​=(1−βn+1​)μ^​n​+βn+1​xn+1​=μ^​n​+βn+1​(xn+1​−μ^​n​)​(1.4)

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

x n + 1 − μ ^ n + 1 = ( x n + 1 − μ ^ n ) + ( μ ^ n − μ ^ n + 1 ) = ( x n + 1 − μ ^ n ) − β n + 1 ( x n + 1 − μ ^ n ) = ( 1 − β n + 1 ) ( x n + 1 − μ ^ n ) (1.5) begin{aligned} 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{aligned} tag{1.5} xn+1​−μ^​n+1​​=(xn+1​−μ^​n​)+(μ^​n​−μ^​n+1​)=(xn+1​−μ^​n​)−βn+1​(xn+1​−μ^​n​)=(1−βn+1​)(xn+1​−μ^​n​)​(1.5)

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

β n ( x n + 1 − μ ^ n + 1 ) 2 = β n ( 1 − β n + 1 ) 2 ( x n + 1 − μ ^ n ) 2 = 1 n n n + 1 n n + 1 ( x n + 1 − μ ^ n ) 2 = β n + 1 ( 1 − β n + 1 ) ( x n + 1 − μ ^ n ) 2 = [ β n + 1 − ( β n + 1 ) 2 ] ( x n + 1 − μ ^ n ) 2 (1.6) begin{aligned} beta _nleft( x_{n+1}-hat{mu}_{n+1} right) ^2&=beta _nleft( 1-beta _{n+1} right) ^2left( 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{aligned} tag{1.6} βn​(xn+1​−μ^​n+1​)2​=βn​(1−βn+1​)2(xn+1​−μ^​n​)2=n1​n+1n​n+1n​(xn+1​−μ^​n​)2=βn+1​(1−βn+1​)(xn+1​−μ^​n​)2=[βn+1​−(βn+1​)2](xn+1​−μ^​n​)2​(1.6)

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

σ ^ n + 1 2 = 1 n ∑ k = 1 n + 1 ( x k − μ ^ n + 1 ) 2 = 1 n ( x n + 1 − μ ^ n + 1 ) 2 + n − 1 n ( 1 n − 1 ∑ k = 1 n [ ( x k − μ ^ n ) + ( μ ^ n − μ ^ n + 1 ) ] 2 ) = β n ( x n + 1 − μ ^ n + 1 ) 2 + ( 1 − β n ) 1 n − 1 ∑ k = 1 n [ ( x k − μ ^ n ) + ( μ ^ n − μ ^ n + 1 ) ] 2 = β n ( x n + 1 − μ ^ n + 1 ) 2 + ( 1 − β n ) 1 n − 1 ∑ k = 1 n ( x k − μ ^ n ) 2 + ( μ ^ n − μ ^ n + 1 ) 2 = [ β n + 1 − ( β n + 1 ) 2 ] ( x n + 1 − μ ^ n ) 2 + ( 1 − β n ) σ ^ n 2 + ( β n + 1 ) 2 ( x n + 1 − μ ^ n ) 2 = ( 1 − β n ) σ ^ n 2 + β n + 1 ( x n + 1 − μ ^ n ) 2 (1.7) begin{aligned} 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 _nleft( 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 _nleft( 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) ^2left( 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{aligned} tag{1.7} σ^n+12​​=n1​k=1∑n+1​(xk​−μ^​n+1​)2=n1​(xn+1​−μ^​n+1​)2+nn−1​(n−11​k=1∑n​[(xk​−μ^​n​)+(μ^​n​−μ^​n+1​)]2)=βn​(xn+1​−μ^​n+1​)2+(1−βn​)n−11​k=1∑n​[(xk​−μ^​n​)+(μ^​n​−μ^​n+1​)]2=βn​(xn+1​−μ^​n+1​)2+(1−βn​)n−11​k=1∑n​(xk​−μ^​n​)2+(μ^​n​−μ^​n+1​)2=[βn+1​−(βn+1​)2](xn+1​−μ^​n​)2+(1−βn​)σ^n2​+(βn+1​)2(xn+1​−μ^​n​)2=(1−βn​)σ^n2​+βn+1​(xn+1​−μ^​n​)2​(1.7)

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

∑ k = 1 n ( x k − μ ^ n ) ( μ ^ n − μ ^ n + 1 ) = ( μ ^ n − μ ^ n + 1 ) ∑ k = 1 n ( x k − μ ^ n ) = ( μ ^ n − μ ^ n + 1 ) ( ∑ k = 1 n x k − ∑ k = 1 n x k ) = 0 (1.8) begin{aligned} 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{aligned}tag{1.8} k=1∑n​(xk​−μ^​n​)(μ^​n​−μ^​n+1​)​=(μ^​n​−μ^​n+1​)k=1∑n​(xk​−μ^​n​)=(μ^​n​−μ^​n+1​)(k=1∑n​xk​−k=1∑n​xk​)=0​(1.8)

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

{ μ ^ n + 1 = μ ^ n + β n + 1 ( x n + 1 − μ ^ n ) σ ^ n + 1 2 = ( 1 − β n ) σ ^ n 2 + β n + 1 ( x n + 1 − μ ^ n ) 2 (1.9) begin{cases} 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{cases}tag{1.9} {μ^​n+1​=μ^​n​+βn+1​(xn+1​−μ^​n​)σ^n+12​=(1−βn​)σ^n2​+βn+1​(xn+1​−μ^​n​)2​(1.9)

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

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      // 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;
}