stmoonar

了解為什麼計算指數加權方差不能正確估計方差。

作者：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$ 的方差！

事實上，如果重複一次上面的 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：