「統計学実践ワークブック」の$8$章「統計的推定の基礎」のバイアスやジャックナイフ推定量に関する直感的な理解が難しいように思われたので、当記事では具体的に正規分布から標本を発生させることでいくつかの値の確認を行いました。ワークブック問$8.3$を主に参考に作成を行いました。
Contents
概要の確認
$\mu^2$の$2$つの推定量$T_1, T_2$の定義
ワークブック問$8.3$では$\mu^2$の推定にあたって、標本$X_1,…,X_n \sim \mathcal{N}(\mu,\sigma^2), i.i.d.$から下記で表す推定量$T_1, T_2$がそれぞれ定義される。
$$
\large
\begin{align}
T_1 &= \bar{X}^{2} \\
T_2 &= \frac{1}{n} \sum_{i=1}^{n} X_{i}^{2}
\end{align}
$$
また、推定量$T_1, T_2$の期待値$E[T_1], E[T_2]$はそれぞれ下記のように計算できる。
$$
\large
\begin{align}
E[T_1] &= E[\bar{X}^2] \\
&= E[\bar{X}]^2 + V[\bar{X}] \\
&= \mu^2 + \frac{\sigma^2}{n} \\
E[T_2] &= E \left[ \frac{1}{n} \sum_{i=1}^{n} X_{i}^{2} \right] \\
&= \frac{1}{n} E \left[ \sum_{i=1}^{n} X_{i}^{2} \right] \\
&= \frac{n}{n} E[X_{i}^{2}] \\
&= E[X_{i}]^{2} + V[X_{i}] \\
&= \mu^{2} + \sigma^{2}
\end{align}
$$
バイアスの定義と推定量$T_1, T_2$のバイアスの計算
母数を$\theta$、推定量を$\hat{\theta}$とおくとき、$\theta$のバイアス$b(\theta)$は下記のように定義される。
$$
\large
\begin{align}
b(\theta) = \hat{\theta} – \theta
\end{align}
$$
上記がわかりにくい場合は「bias = parameter – estimator」のように考えると良い。ここで前項のように定めた推定量$T_1,T_2$のバイアスを$b_{1}(\mu,\sigma^2),b_{2}(\mu,\sigma^2)$のようにおくと、それぞれ下記のように計算できる。
$$
\large
\begin{align}
b_{1}(\mu,\sigma^2) &= E[T_1] – \mu^2 \\
&= \mu^2 + \frac{\sigma^2}{n} – \mu^2 \\
&= \frac{\sigma^2}{n} \\
b_{2}(\mu,\sigma^2) &= E[T_2] – \mu^2 \\
&= \mu^2 + \sigma^2 – \mu^2 \\
&= \sigma^2
\end{align}
$$
ここでの導出は「統計学実践ワークブック」の問$8.3$の$[1]$に対応する。また、式を確認することで$T_1$のバイアスは$n$が大きいときには$0$に収束する一方で、$T_2$のバイアスは一定であることが確認できる。
ジャックナイフ推定量
詳しくは「統計学実践ワークブック」の$8$章まとめの「ジャックナイフ推定量」で取り扱ったように、ジャックナイフ推定量はバイアスの推定量を元に下記のように定められる。
$$
\large
\begin{align}
\hat{\theta}_{(\cdot)} &= \frac{1}{n} \sum_{i=1}^{n} \hat{\theta}_{(i)} \\
\hat{\mathrm{bias}} &= (n-1)(\hat{\theta}_{(\cdot)}-\hat{\theta}) \\
\hat{\theta}_{jack} &= \hat{\theta} – \hat{\mathrm{bias}}
\end{align}
$$
$\hat{\theta}_{(\cdot)}$の計算手順が少々複雑であるので注意が必要である。推定量$T_1$にあてはめる場合は、$i$番目の標本$X_i$を除いて計算を行なった標本平均を$\bar{X}_{(i)}$とおき、$\hat{\theta}_{(i)}=\bar{X}_{(i)}^2, \hat{\theta}=\bar{X}^2$を元に計算を行えばよい。
Pythonを用いた具体的な値の確認
バイアスの収束の確認
「バイアスの定義と推定量$T_1, T_2$のバイアスの計算」で取り扱った内容に関して、以下Pythonを用いて$\mathcal{N}(10,1^2)$を元に生成した乱数に関して確認を行う。
import numpy as np
from scipy import stats
np.random.seed(0)
mu, sigma, num_trial = 10., 1., 50
num_sample = np.array([5,10,20,100,500,1000,10000,100000])
for i in range(num_sample.shape[0]):
T1_sum, T2_sum = 0., 0.
for j in range(num_trial):
x = stats.norm.rvs(mu,sigma,size=num_sample[i])
T1_sum += np.mean(x)**2
T2_sum += np.mean(x**2)
T1, T2 = T1_sum/num_trial, T2_sum/num_trial
print("---")
print("Num_Sample: {:.0f}, T1-mu^2: {:.2f}".format(num_sample[i],T1-mu**2))
print("Num_Sample: {:.0f}, T2-mu^2: {:.2f}".format(num_sample[i],T2-mu**2))
・実行結果
Num_Sample: 5, T1-mu^2: 0.87
Num_Sample: 5, T2-mu^2: 1.63
---
Num_Sample: 10, T1-mu^2: -2.34
Num_Sample: 10, T2-mu^2: -1.40
---
Num_Sample: 20, T1-mu^2: 0.61
Num_Sample: 20, T2-mu^2: 1.50
---
Num_Sample: 100, T1-mu^2: -0.29
Num_Sample: 100, T2-mu^2: 0.69
---
Num_Sample: 500, T1-mu^2: -0.02
Num_Sample: 500, T2-mu^2: 0.97
---
Num_Sample: 1000, T1-mu^2: -0.03
Num_Sample: 1000, T2-mu^2: 0.96
---
Num_Sample: 10000, T1-mu^2: 0.07
Num_Sample: 10000, T2-mu^2: 1.07
---
Num_Sample: 100000, T1-mu^2: 0.00
Num_Sample: 100000, T2-mu^2: 1.00
上記の実行結果より、概ね$500$サンプル以上を用いて計算を行なった場合、$\displaystyle b_{1}(\mu,\sigma^2)=\frac{\sigma^2}{n}$が$0$に近づき、$b_{2}(\mu,\sigma^2)=\sigma^2$が$0$に近づくことが読み取れる。
ジャックナイフ法を用いた$T_1$の補正
以下では、「ジャックナイフ推定量」で取り扱った内容を元にジャックナイフ法を用いて$T_1$の補正を行い、具体的な値の変化の確認を行う。
import numpy as np
from scipy import stats
np.random.seed(10)
mu, sigma, num_trial = 10., 1., 5
num_sample = np.array([3,5,10,15,20,50,100])
for i in range(num_sample.shape[0]):
T1_sum, T1_jack_sum = 0., 0.
for j in range(num_trial):
x = stats.norm.rvs(mu,sigma,size=num_sample[i])
x_sub = np.zeros(num_sample[i])
for k in range(num_sample[i]):
a = 0
for l in range(num_sample[i]):
if k != l:
a += x[l]
x_sub[k] = (a/(x.shape[0]-1))**2
x_sub_mean = np.mean(x_sub)
bias = (x.shape[0]-1)*(x_sub_mean - np.mean(x)**2)
T1_sum += np.mean(x)**2
T1_jack_sum += np.mean(x)**2-bias
T1, T1_jack = T1_sum/num_trial, T1_jack_sum/num_trial
print("---")
print("Num_Sample: {:.0f}, T1-mu^2: {:.2f}".format(num_sample[i],T1-mu**2))
print("Num_Sample: {:.0f}, T1_jack-mu^2: {:.2f}".format(num_sample[i],T1_jack-mu**2))
・実行結果
Num_Sample: 3, T1-mu^2: 3.43
Num_Sample: 3, T1_jack-mu^2: 3.14
---
Num_Sample: 5, T1-mu^2: 3.58
Num_Sample: 5, T1_jack-mu^2: 3.32
---
Num_Sample: 10, T1-mu^2: 0.23
Num_Sample: 10, T1_jack-mu^2: 0.13
---
Num_Sample: 15, T1-mu^2: 1.74
Num_Sample: 15, T1_jack-mu^2: 1.66
---
Num_Sample: 20, T1-mu^2: -0.97
Num_Sample: 20, T1_jack-mu^2: -1.01
---
Num_Sample: 50, T1-mu^2: 0.79
Num_Sample: 50, T1_jack-mu^2: 0.77
---
Num_Sample: 100, T1-mu^2: -1.16
Num_Sample: 100, T1_jack-mu^2: -1.16
実行結果より、ジャックナイフ推定量を用いることで推定された値が誤差が小さくなるように補正がかけられていることが確認できる。また、ここでバイアスは概ね正の値を取ることも合わせて確認できた。
[…] ・考察ここで行なった結果が妥当かどうか確認するにあたって、下記で正規分布からの標本を元にPythonを用いて結果の検証を行なった。https://www.hello-statisticians.com/explain-terms-cat/bias1.html […]