了解为什么计算指数加权方差不能正确估计方差。
作者: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)。
估计量的偏差是指估计量的预期值与被估计参数(本例中为方差)的真实值之间的差值。有偏差的估计值与真实值之差不为零,而无偏差的估计值与真实值之差为零。
让我们试着测量方差,看看会发生什么。随机变量 $$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 个随机点的平均值。
事实上,如果重复一次上面的 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 指数加权方差所需的所有工具。
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: