区間推定と仮説検定は推測統計を学ぶ上で主要なトピックです。一方で、教科書や参考書の解説は抽象的で直感的な理解が難しい場合も多いです。そこで当記事では$95$％区間が本当に$95$％で成立するかに関してSciPyを用いて生成したサンプルを元に検証実験を行いました。
・コラム
https://www.hello-statisticians.com/column
基本方針の確認
サンプルの生成
サンプルの生成はscipy.statsを用いて行うことができます。たとえば標準正規分布$\mathcal{N}(0,1)$に従う$n=20$のサンプルは下記のように生成できます。
import numpy as np
from scipy import stats
np.random.seed(10)
x = stats.norm.rvs(size=20)
print(x)・実行結果
[ 1.3315865 0.71527897 -1.54540029 -0.00838385 0.62133597 -0.72008556
0.26551159 0.10854853 0.00429143 -0.17460021 0.43302619 1.20303737
-0.96506567 1.02827408 0.22863013 0.44513761 -1.13660221 0.13513688
1.484537 -1.07980489]上記のプログラムではnp.random.seed(10)を実行することで乱数のseedを固定しました。seedが変わると結果の再現ができなくなるので、検証実験などを行う場合は必須になります。
実験の方針
区間推定によって得られる$95$％区間は「$95$％の確率で推定した区間に母数(パラメータ)が含まれる」ことを意味します。もう少し直感的に考えると「$20$回推定を行った場合に$19$回は正しい区間が得られ、$1$回は正しくない区間が得られる」と解釈することができます。
よって、次節における実験では正規分布$\mathcal{N}(1,3^2)$に基づいて$n=20$のサンプルを元に区間推定を$300$回行った際に、正しい区間が何回得られるかに関して検証を行います。
検証実験
母分散未知の場合の母平均の推定
$$
\large
\begin{align}
t_{\alpha=0.975}(n-1) \leq & \frac{\sqrt{n}(\overline{X}-\mu)}{s} \leq t_{\alpha=0.025}(n-1) \\
-\frac{s}{\sqrt{n}} t_{\alpha=0.025}(n-1) \leq & \mu-\overline{X} \leq -\frac{s}{\sqrt{n}} t_{\alpha=0.975}(n-1) \\
\overline{X} – \frac{s}{\sqrt{n}} t_{\alpha=0.025}(n-1) \leq & \mu \leq \overline{X} + \frac{s}{\sqrt{n}} t_{\alpha=0.025}(n-1)
\end{align}
$$
上側$2.5$％点を$t_{\alpha=0.025}$のように表すとき、自由度$n-1$の$95$％区間は上記のように表すことができます。上記が$95$％の確率で成立するかを以下を実行することで検証を行います。
import numpy as np
from scipy import stats
np.random.seed(15)
alpha = 0.025
for i in range(300):
x = stats.norm.rvs(loc=1, scale=3, size=20)
s2 = np.sum((x-np.mean(x))**2)/(20.-1.)
L_mu, U_mu = np.mean(x)-np.sqrt(s2/20.)*stats.t.ppf(1-0.025,20-1), np.mean(x)+np.sqrt(s2/20.)*stats.t.ppf(1-0.025,20-1)
if 1 < L_mu or U_mu < 1:
print("Experiment.{}, L_mu: {:.3f}, U_mu: {:.3f}".format(i+1, L_mu, U_mu))・実行結果
Experiment.81, L_mu: -0.974, U_mu: 0.969
Experiment.91, L_mu: 1.390, U_mu: 3.136
Experiment.114, L_mu: 1.035, U_mu: 3.295
Experiment.122, L_mu: 1.390, U_mu: 4.255
Experiment.138, L_mu: -1.843, U_mu: 0.548
Experiment.158, L_mu: 1.007, U_mu: 2.841
Experiment.190, L_mu: -1.916, U_mu: 0.834
Experiment.212, L_mu: -3.025, U_mu: 0.544
Experiment.227, L_mu: 1.849, U_mu: 4.206
Experiment.234, L_mu: 1.390, U_mu: 4.532
Experiment.235, L_mu: 2.005, U_mu: 4.980
Experiment.236, L_mu: 1.205, U_mu: 4.013
Experiment.239, L_mu: 1.292, U_mu: 3.212
Experiment.292, L_mu: 1.039, U_mu: 3.742
Experiment.298, L_mu: -1.644, U_mu: 0.980上記より$300$回中$15$回の区間推定が正しくない結果を出力したことが確認できます。よって、乱数のseedに注意が必要ですが、母平均の区間推定に関して$95$％区間が概ね$95$％で成立するであろうことが上記から考察できます。
母分散の推定
$$
\large
\begin{align}
\chi^{2}_{\alpha=0.975}(n-1) \leq & \frac{(n-1)s^2}{\sigma^2} \leq \chi^{2}_{\alpha=0.025}(n-1) \\
\frac{1}{\chi^{2}_{\alpha=0.025}(n-1)} \leq & \frac{\sigma^2}{(n-1)s^2} \leq \frac{1}{\chi^{2}_{\alpha=0.975}(n-1)} \\
\frac{(n-1)s^2}{\chi^{2}_{\alpha=0.025}(n-1)} \leq & \sigma^2 \leq \frac{(n-1)s^2}{\chi^{2}_{\alpha=0.975}(n-1)}
\end{align}
$$
上側$2.5$％点を$\chi^{2}_{\alpha=0.025}$のように表すとき、自由度$n-1$の$95$％区間は上記のように表すことができます。上記が$95$％の確率で成立するかを以下を実行することで検証を行います。
import numpy as np
from scipy import stats
np.random.seed(10)
alpha = 0.025
for i in range(300):
x = stats.norm.rvs(loc=1, scale=3, size=20)
s2 = np.sum((x-np.mean(x))**2)/(20.-1.)
L_sigma2, U_sigma2 = (20.-1.)*s2/stats.chi2.ppf(1-alpha,20-1), (20.-1.)*s2/stats.chi2.ppf(alpha,20-1)
if 3**2 < L_sigma2 or U_sigma2 < 3**2:
print("Experiment.{}, L_sigma2: {:.3f}, U_sigma2: {:.3f}".format(i+1, L_sigma2, U_sigma2))・実行結果
Experiment.12, L_sigma2: 2.437, U_sigma2: 8.989
Experiment.20, L_sigma2: 1.820, U_sigma2: 6.712
Experiment.28, L_sigma2: 1.273, U_sigma2: 4.697
Experiment.30, L_sigma2: 9.716, U_sigma2: 35.840
Experiment.62, L_sigma2: 9.903, U_sigma2: 36.529
Experiment.66, L_sigma2: 2.363, U_sigma2: 8.717
Experiment.83, L_sigma2: 2.007, U_sigma2: 7.403
Experiment.131, L_sigma2: 10.366, U_sigma2: 38.235
Experiment.138, L_sigma2: 2.439, U_sigma2: 8.998
Experiment.171, L_sigma2: 9.448, U_sigma2: 34.848
Experiment.201, L_sigma2: 9.716, U_sigma2: 35.837
Experiment.277, L_sigma2: 0.971, U_sigma2: 3.580
Experiment.281, L_sigma2: 10.045, U_sigma2: 37.053
Experiment.285, L_sigma2: 1.894, U_sigma2: 6.984
Experiment.287, L_sigma2: 2.186, U_sigma2: 8.062上記より$300$回中$15$回の区間推定が正しくない結果を出力したことが確認できます。よって、乱数のseedに注意が必要ですが、母分散の区間推定に関して$95$％区間が概ね$95$％で成立するであろうことが上記から考察できます。