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）です。

バイアスとは何ですか？#

推定量のバイアスは、推定量の期待値と推定されるパラメータ（この場合は分散）の真の値との差を指します。有偏推定値は真の値との差がゼロではなく、無偏推定値は真の値との差がゼロです。

分散を測定してみましょう。確率変数 $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  # サンプル数

# これにより、各行が1つのシミュレーションで、列がサンプル値の配列が得られます。
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 は、中心点（すなわち平均）が 0 の 100 個のランダムな点の例を示しています。その中から 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}$

我々の推定値は $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]$ の正確な値を計算します。展開することでこれを実現できます。次のように始めます：

$\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é の定理により、これは次のように等しいことがわかります：

$\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

# フェイクデータを作成
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

    # 重み
    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 が得られます：