复现Pandas指数加权方差

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

作者：Adrian Letchford 译者：stmoonar

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

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

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

1

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

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

什么是偏差（bias）？

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

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

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

然后用样本的平均值代替 $\mu$：

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

 1
 2
 3
 4
 5
 6
 7
 8
 9
10
11
12
13

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。那么我们的平均估计值是多少呢？计算估计值的平均值：

1

estimates.mean()

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

偏差从哪里来？

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

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

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

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

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

1

avg = 0

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

1

estimates.mean()

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

量化偏差

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

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

开始。

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

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

现在，请注意：

 1
 2
 3
 4
 5
 6
 7
 8
 9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38

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