正規分布からの標本に基づく母平均の二乗$\mu^2$の推定におけるバイアスとジャックナイフ推定量

「統計学実践ワークブック」の$8$章「統計的推定の基礎」のバイアスやジャックナイフ推定量に関する直感的な理解が難しいように思われたので、当記事では具体的に正規分布から標本を発生させることでいくつかの値の確認を行いました。ワークブック問$8.3$を主に参考に作成を行いました。

日本統計学会公式認定統計検定準1級対応統計学実践ワークブック

3,080円(04/20 02:24時点)

Amazon

Contents

1 概要の確認
2 Pythonを用いた具体的な値の確認

概要の確認

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