复现 Pandas 指数加权方差

了解为什么计算指数加权方差不能正确估计方差。

作者：Adrian Letchford 译者：stmoonar

你很可能熟悉用指数加权法计算平均数的方法。指数加权平均法的原理是对较新的信息赋予较高的权重。指数加权平均法的权重如下：

$\displaystyle{w_t = (1 - \alpha)^t}$

对于 $t \in [0, \dots, T]$ ，序列 $X_t$ 的指数加权平均值是这样的：

$\displaystyle{\bar{X}_T = \frac{1}{\sum_t w_t} \sum_t w_t X_t}$

在 Pandas 中，您可以通过以下面的方法轻松计算指数加权移动平均线：

df.ewm(alpha=0.1).mean()

如果你自己计算指数加权平均值，你会发现它与 Pandas 给出的结果一致。但是，我们即将看到，如果我们尝试用方差来计算，我们将得到一个很差的估计值。这就是所谓的估计偏差（estimation bias）。

什么是偏差（bias）？#

估计量的偏差是指估计量的预期值与被估计参数（本例中为方差）的真实值之间的差值。有偏差的估计值与真实值之差不为零，而无偏差的估计值与真实值之差为零。

让我们试着测量方差，看看会发生什么。随机变量 $X$ 的方差是：

$\displaystyle{\sigma^2 = E[(X - \mu)^2]}$

如果我们有 $X$ 的 $n$ 值样本，我们可以尝试用样本的平均值代替期望值 $E[\cdot]$ 来估计方差：

$\displaystyle{\frac{1}{n} \sum_i \left(X_i - \mu \right)^2}$

然后用样本的平均值代替 $\mu$ ：
$\displaystyle{\begin{aligned} \bar{X} &= \frac{1}{n}\sum_i X_i \\ \hat{\sigma}^2 &= \frac{1}{n} \sum_i \left(X_i - \bar{X} \right)^2 \end{aligned}}$

我们可以用 Python 写一个简单的模拟，看看我们的估计量得到的 $\hat{\sigma}^2$ 有多好。

from scipy.stats import norm

num_simulations = 10_000_000
n = 5  # 样本数

# 这样就得到了一个数组，其中每一行都是一个模拟，而列则是样本值。
simulations = norm.rvs(loc=0, scale=1, size=(num_simulations, n))

# 由此得出每个模拟的估计平均值。
avg = simulations.mean(axis=1).reshape(-1, 1)

# 使用我们的估计器来估计每个模拟的方差。
estimates = ((simulations - avg)**2).sum(1) / n

译者注：

代码使用了 norm.rvs 方法从正态分布中随机生成样本值，loc=0 表示正态分布的均值为 0，scale=1 表示标准差为 1，size=(num_simulations, n) 表示生成的数组将有 num_simulations 行，每行有 n 个样本值。

现在我们有 1000 万个方差估计值，而真实方差为 1。那么我们的平均估计值是多少呢？计算估计值的平均值：

estimates.mean()

得出约 0.8！我们的估计值偏差了 0.2。这就是偏差！

偏差从哪里来？#

让我们回到方差的定义以及如何将其转化为样本估计值。为了计算我们的估计值，我们将 $\mu$ 换成了样本平均值：

$\displaystyle{\frac{1}{n} \sum_i \left(X_i - \mu \right)^2 \quad\Rightarrow\quad \frac{1}{n} \sum_i \left(X_i - \bar{X} \right)^2}$

这就是偏差产生的原因。小样本的平均值（ $\bar{X}$ ）会比总体的平均值（ $\mu$ ）更接近这些样本，通过下面的例子可以更加直观。

图 1 显示了一个有 100 个随机点的例子，这些点的中心点（即平均值）为 0，其中有 5 个点被随机选中，它们的平均值用黑叉表示。这 5 个样本的平均值就是最接近这 5 个样本的点。根据定义，样本的平均数会比总体平均数更接近样本。因此：

$\displaystyle{\frac{1}{n} \sum_i \left(X_i - \bar{X} \right)^2 \quad\lt\quad \frac{1}{n} \sum_i \left(X_i - \mu \right)^2}$

我们的估计值会低于（underestimate ） $X$ 的方差！

图 1：偏差示意图。这是一幅平均值为（0，0）的 100 个随机点的图。图中随机选取了 5 个点并突出显示。黑叉表示这 5 个随机点的平均值。

图 1：偏差示意图。 这是一幅平均值为（0，0）的 100 个随机点的图。图中随机选取了 5 个点并突出显示。黑叉表示这 5 个随机点的平均值。

事实上，如果重复一次上面的 Python 模拟，但是使用 0（总体平均值）代替样本平均值：

avg = 0

那么得到的样本方差的平均值就是 1：

estimates.mean()

如果知道总体均值，我们就可以从一组样本中得到方差的无偏估计值。但实际上，我们并不知道总体均值。幸运的是，我们可以量化偏差并加以纠正。

量化偏差#

到目前为止，我们已经看到 $\hat{\sigma}^2$ 是对总体方差的一个有偏差的估计。我们通过模拟多个样本得到 $\hat{\sigma}^2$ 并取平均值发现了这一点。模拟结果表明：

$\displaystyle{E[\hat{\sigma}^2] \ne \sigma^2}$

我们现在要摆脱模拟，计算 $E[\hat{\sigma}^2]$ 的精确值。我们可以通过展开式（原文：expanding it out）来实现。我们从：

$\displaystyle{E[\hat{\sigma}^2] = E \left[ \frac{1}{n} \sum_i \left(X_i - \bar{X} \right)^2 \right]}$

开始。

我们可以让 $\bar{X} = \mu - (\mu - \bar{X})$ ，这意味着我们可以展开为

$\displaystyle{E[\hat{\sigma}^2] = E \left[ \frac{1}{n} \sum_i \left((X_i - \mu)- (\bar{X} - \mu) \right)^2 \right]}$

通过一些代数知识，我们可以展开平方式：

$\displaystyle{\begin{aligned} E[\hat{\sigma}^2] &= E \left[ \frac{1}{n} \sum_i \left((X_i - \mu)^2 - 2(\bar{X} - \mu)(X_i - \mu) + (\bar{X} - \mu)^2\right) \right] \\ &= E \left[ \frac{1}{n} \sum_i (X_i - \mu)^2 - 2(\bar{X} - \mu) \frac{1}{n} \sum_i(X_i - \mu) + \frac{1}{n} \sum_i(\bar{X} - \mu)^2 \right] \\ &= E \left[ \frac{1}{n} \sum_i (X_i -\mu)^2 - 2(\bar{X} - \mu) \frac{1}{n} \sum_i(X_i - \mu) + (\bar{X} -\mu)^2 \right] \end{aligned}}$

现在，请注意：

$\displaystyle{\frac{1}{n} \sum_i(X_i - \mu) = \frac{1}{n} \sum_i X_i - \frac{1}{n} \sum_i \mu = \frac{1}{n} \sum_i X_i - \mu = \bar{X} - \mu}$

这意味着：

$\displaystyle{\begin{aligned} E[\hat{\sigma}^2] &= E \left[ \frac{1}{n} \sum_i (X_i - \mu)^2 - 2(\bar{X} - \mu)^2 + (\bar{X} - \mu)^2 \right] \\ &= E \left[ \frac{1}{n} \sum_i (X_i - \mu)^2 - (\bar{X} - \mu)^2 \right] \\ &= E \left[ \frac{1}{n} \sum_i (X_i - \mu)^2 \right] - E \left[ (\bar{X} - \mu)^2 \right] \end{aligned}}$

这里的妙处是：

$\displaystyle{E \left[ \frac{1}{n} \sum_i (X_i - \mu)^2 \right] = \sigma^2}$

意味着：

$\displaystyle{\begin{aligned} E[\hat{\sigma}^2] &= \sigma^2 - E \left[ (\bar{X} - \mu)^2 \right] \end{aligned}}$

项 $E \left[ (\bar{X} - \mu)^2 \right]$ 是样本平均数的方差。根据Bienaymé’s identity ，我们知道它等于

$\displaystyle{E \left[ (\bar{X} - \mu)^2 \right] = \frac{1}{n}\sigma^2}$

这让我们得到：

$\displaystyle{E[\hat{\sigma}^2] = (1 - \frac{1}{n}) \sigma^2 = (1 - \frac{1}{5}) \times 1 = 0.8}$

回想一下我们的 Python 模拟；样本数为 $n=5$ ，真实方差为 $\sigma^2 = 1$ ，估计方差为 $\hat{\sigma}^2 = 0.8$ 。如果我们将 $n$ 和 $\sigma^2$ 插入上述公式，就会得到有偏差的答案：

$\displaystyle{E[\hat{\sigma}^2] = (1 - \frac{1}{n}) \sigma^2 = (1 - \frac{1}{5}) \times 1 = 0.8}$

无偏估计量#

现在我们知道了 $E[\hat{\sigma}^2]$ 的精确值，就可以想办法修正 $\hat{\sigma}^2$ ，使其成为 $\sigma^2$ 的无偏估计值。

修正项为：

$\displaystyle{\frac{n}{n-1}}$

通过演算，我们可以看到：

$\displaystyle{\begin{aligned} E[\frac{n}{n-1} \hat{\sigma}^2] &=\frac{n}{n-1} E[\hat{\sigma}^2] \\ &= \frac{n}{n-1}(1 - \frac{1}{n}) \sigma^2 \\ &= \frac{n(1 - \frac{1}{n})}{n-1} \sigma^2 \\ &= \frac{n - 1}{n-1} \sigma^2 \\ &= \sigma^2 \end{aligned}}$

因此，从一组样本中得到的 $\sigma^2$ 的无偏估计值是：

$\displaystyle{\begin{aligned} \frac{n}{n-1} \hat{\sigma}^2 &= \frac{n}{n-1} \frac{1}{n} \sum_i \left(X_i - \bar{X} \right)^2 \\ &= \frac{1}{n-1} \sum_i \left(X_i - \bar{X} \right)^2 \end{aligned}}$

无偏加权估计量#

现在，我们将上述内容扩展到样本加权的情况。加权样本平均值为：

$\displaystyle{\bar{X} = \frac{1}{\sum_i w_i} \sum_i w_i X_i}$

加权方差为：

$\hat{\sigma}^2 = \frac{1}{\sum_i w_i} \sum_i w_i\left(X_i - \bar{X} \right)^2$

按照与之前完全相同的展开过程，我们得出：

$\displaystyle{\begin{aligned} E[\hat{\sigma}^2] &= \sigma^2 - E \left[ (\bar{X} - \mu)^2 \right] \end{aligned}}$

平均值的方差为：

$\displaystyle{\begin{aligned} E \left[ (\bar{X} - \mu)^2 \right] &= \text{var}(\bar{X}) \\ &= \text{var}\left(\frac{1}{\sum w_i} \sum w_i X_i \right) \\ &= \frac{1}{(\sum w_i)^2} \sum \text{var} (w_i X_i) \\ &= \frac{1}{(\sum w_i)^2} \sum w_i^2 \text{var} (X_i) \\ &= \frac{\sum w_i^2}{(\sum w_i)^2} \sigma^2 \end{aligned}}$

$var$ 是计算方差。

这样我们就得到：

$\displaystyle{\begin{aligned} E[\hat{\sigma}^2] &= \sigma^2 - \frac{\sum w_i^2}{(\sum w_i)^2} \\ \sigma^2 &= \left(1 - \frac{\sum w_i^2}{(\sum w_i)^2} \right)\sigma^2 \end{aligned}}$

偏差修正项为：

$\displaystyle{b = \frac{(\sum w_i)^2}{(\sum w_i)^2 - \sum w_i^2}}$

这意味着方差的无偏加权估计值为：

$\displaystyle{b \hat{\sigma}^2}$

复现 Pandas 指数加权方差#

现在，我们拥有了复现 Pandas 指数加权方差所需的所有工具。

import numpy as np
import pandas as pd

N = 1000

# 创建一些fake data
df = pd.DataFrame()
df['data'] = np.random.randn(N)

# 为 EWM 设置半衰期，并将其转换为 alpha 值进行计算。
halflife = 10
a = 1 - np.exp(-np.log(2)/halflife)  # alpha

# 这是 Pandas 的 ewm
df['var_pandas'] = df.ewm(alpha=a).var()

# 初始化变量
varcalc = np.zeros(len(df))

# 计算指数移动方差
for i in range(0, N):

    x = df['data'].iloc[0:i+1].values

    # Weights
    n = len(x)
    w = (1-a)**np.arange(n-1, -1, -1) # 顺序相反，以便与序列顺序一致

    # 计算指数移动平均线
    ewma = np.sum(w * x) / np.sum(w)

    # 计算偏差修正项
    bias = np.sum(w)**2 / (np.sum(w)**2 - np.sum(w**2))

    # 计算带偏差修正的指数移动方差
    varcalc[i] = bias * np.sum(w * (x - ewma)**2) / np.sum(w)

df['var_calc'] = varcalc

这样我们就得到了一个下面这样的 DataFrame：