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

文章目录

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

前言

对于一个序列而言，求均值和方差根据定义式是不难的，其时空复杂度均为 O ( N ) \mathcal{O}\left(N\right) O(N) 。但有的时候，我们的样本是一个一个给的，此时新来了一个样本，我们总不可能把原来的样本都捞出来再算一次均值、方差吧，那样时空复杂度都是 O ( N ) \mathcal{O}\left(N\right) 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(1\right) 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.4\right) (1.4) 稍微变形，引入式 ( 1.5 ) \left(1.5\right) (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.6\right) (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 _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{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.7\right) (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 _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{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.7\right) (1.7) 的推导第三行到第四行的等号是因为完全平方展开的交叉项为 0 0 0 ，具体推导可见式 ( 1.8 ) \left(1.8\right) (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)

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