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  # サンプル数

# これにより、各行がシミュレーションで、列がサンプル値の配列が得られます。
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 が得られます：