とても有名な「統計学入門」をpythonで実装しながら読んでいきます。
本をしっかり読むことと実装の勉強が目的ですので説明が不足していたりします。
また、実装方法に関してはベストなものではないと思いますのでご了承ください。

前回
正規分布からの標本

推定

母集団の母数は実際の問題では未知であり、これらを $X_1,X_2,\cdots,X_n$ から定める必要がある。
これを母数の推定と呼ぶ。
標本平均や標本分散のように、母数を推定するために標本から求められた統計量を一般に推定量と呼ぶ。
推定しようとするパラメータを$\theta$で表すものとする。推定量は、$\hat{\theta}$のように表す。
たとえば、$\theta$が母平均のとき、$\theta$を推定するために、標本平均を

\hat{\theta}=(X_1+\cdots+X_n)\ /\ n

のように使う。
一般に$\hat{\theta}$は$X_1,\cdots,X_n$の関数であるが、それを強調する場合は$\hat{\theta}(X_1,\cdots,X_n)$と書くこともある。
また、複数個の母数を同時に考える場合は、これらを$\theta_1,\cdots,\theta_k$、その推定量を$\hat{\theta}_1,\cdots,\hat{\theta}_k$と表す。

点推定と区間推定

母集団の未知の母数$\theta$を推定する場合、それをある一つの値$\hat{\theta}$を指定する方法を点推定と呼ぶ。
$\hat{\theta}$は$X_1,\cdots,X_n$の関数である。この関数が推定量であるが、確率変数であるから現実の$\hat{\theta}$の値は$\theta$に一致せず、実際の推定にはなにがしかの誤差を伴う。したがって、誤差の評価に行わなければならないが、これには推定量の標本分布を考えた確率的な取り扱いが必要である。
区間推定は、真のパラメータに値が入る確率がある値$1-\alpha$以上と保証される区間$[L,U]$を求めるもので、最初からある程度の誤りがあることを認めた推定法といえる。たとえば、

P(L\leq \mu \leq U) \geq 1-\alpha

となる、$L,U$を求めるものである。
$L,U$は$X_1,\cdots,X_n$の関数、つまり統計量で、左辺の確率が$1-\alpha$以上になるように標本分布から求められる

点推定の考え方とその手順

たとえば、標本平均$\bar{X}$は母平均$\mu$の推定量、標本分散$s^2$は母分散$\sigma^2$の推定量である。

E(\bar{X})=\mu,\hspace{5mm}E[s^2]=\sigma^2

であるから、平均的に$\bar{X}$は$\mu$に、$s^2$は$\sigma^2$に合致する推定量である。推定量は$X_1,\cdots,X_n$の関数で確率変数である。
具体的に$n$個の観測値が与えられた場合、推定量を実際の数字として計算することができる。これを推定値と呼ぶ。

点推定の手順

モーメント法

母集団が正規分布$N(\mu,\sigma^2)$のとき、母集団の1次、2次モーメント$\mu_1,\mu_2$は$\sigma^2=\mu_2-\mu_1^2$だから、

\mu_1=\mu_2,\hspace{5mm}\mu_2=\sigma^2+\mu^2

である。一方、標本から、

\hat{\mu}_1=\sum X_i\ /\ n,\hspace{5mm}\hat{\mu}_2=\sum X_i^2\ /\ n

で、推定することを考える。母モーメント=標本モーメントと等しく置くと、

\mu_1=\hat{\mu}_1,\hspace{5mm}\mu_2=\hat{\mu}_2

から、母平均、母分散の$\mu,\sigma^2$の推定量

\hat{\mu}=\sum X_i\ /\ n\hspace{5mm}\hat{\sigma}^2=\sum(X_i-\bar{X})^2\ /\ n

が得られる。
このようにモーメント$\mu_1,\mu_2,\cdots$の推定によって、母数の推定を行う手続きをモーメント法という。

最尤法

1をとる確率が$p$、0をとる確率が$1-p$のベルヌーイ分布$Bi(1,p)$が母集団分布の場合を考える。
推定すべき未知の母数は$p$であり、0から1までの間の値をとる可能性があるが、母数が取り得る値の集合を母数空間と呼び、$\Theta$で表す。
たとえば、$X_1=1,X_2=1,X_3=1,X_4=1,X_5=0$の標本が得られたとする。
最尤原理といわれる「現実の標本は確率最大のものが実現した」という仮定を採用する。
この標本が得られる確率は、

L(p)=p^4(1-p)

である。
$L(p)$は、母数空間$\Theta$での$p$のいろいろな値におけるもっともらしさを表す関数とみなすことができ、このようにもっともらしさを尤度、その関数を尤度関数と呼ぶ。
最尤法は尤度関数を母数空間$\Theta$で最大にするものを推定値や推定量とするもので、尤度関数を最大にする値が最尤推定値、関数としては最尤推定量である。
今回の例では、

\frac{dL(p)}{dp}=p^3(4-5p)

なので、$\hat{p}=0.8$が最尤推定値である。

一般に、大きさ$n$の標本において、未知の母数を一般に$\theta$とした場合、尤度関数は$X_1,\cdots,X_n$の同時確率分布を$\theta$の関数とみなしたもので、$X_i$の確率分布を$f(x_i,\theta)$とすると、尤度関数は$n$個の確率関数の積

L(\theta)=f(x_1,\theta)・f(x_2,\theta)・\cdots・f(x_n,\theta)

を$\theta$の関数と見なしたものである。
また、母集団分布が複数の母数$\theta_1,\cdots,\theta_k$を含む場合は、尤度関数は、

L(\theta_1,\cdots,\theta_k)=\prod_{i=1}^nf(x_i;\theta_1,\cdots,\theta_k)

となる。
対数をとって和の形にした対数尤度

\log L(p)=\sum_{i=1}^nf(X_i;\theta)

を考えるのが普通である。
最尤推定量$\hat{\theta}$は、$\log L(\theta)$を母数空間$\Theta$において最大にする推定量であり、通常

\frac{\partial L(\theta)}{\partial \theta}=0

の解を考えることによって求める。
ベルヌーイ分布では、

\log L(p)=r\log p+(n-r)\log (1-p)

これを微分して

\frac{\partial L(p)}{\partial p}=\frac{r}{p}-\frac{n-r}{1-p}=0

とおくと、$p$の最尤推定量は1の相対頻度$\hat{p}=r\ /\ n$となる。

点推定の基準

推定量は一つとは限らず複数考えることができる。
すべてが推定量として適当であるとは限らず、実用上意味があるためには、標本分布は真の母数$\theta$の周辺に集中していることを示すいくつかの基準を満足する必要がある。

基準1 不偏性

$\hat{\theta}$は$\theta$の推定量であるから、その分布は$\theta$の回りに散布していなければならない。
この一つの基準が不偏性である。
不偏性とは、推定量の期待値をとった場合、真の母数の値となることであり、

E[\hat{\theta}]=\theta

となる性質である。これを満足する推定量を不偏推定量と呼ぶ。
いま、母平均を$\mu$、母分散を$\sigma^2$とすると、大きさ$n$の無作為標本$X_1,\cdots,X_n$に対し、$E(X_i)=\mu,V(X_i)=\sigma^2$である。
標本平均

\bar{X}=(X_1+\cdots+X_n)\ /\ n

に対しては、

E(\bar{X})=\mu

だから、標本平均はつねに母平均の不偏推定量である。
母分散$\sigma^2$の推定のために、

S^2=\frac{1}{n}\sum_{i=1}^n(X_i-\bar{X})^2

を使うと、計算の結果

E(S^2)=\frac{n-1}{n}\sigma^2

となり、平均的に過小推定することになり、$\sigma^2$に対して不偏ではない。

s^2=\frac{1}{n-1}\sum_{i=1}^n(X_i-\bar{X})^2

を用いれば、

E(s^2)=\sigma^2

となって、母分散の不偏推定量$\sigma^2$を得る。

基準２一致性

標本の大きさ$n$が大きくなるに従い、$\hat{\theta}_n$が、真の母数の値$\theta$に近づく性質である。
数学的には、すべての$\varepsilon>0$に対して、$n→\infty$のとき$P(|\hat{\theta}_n-\theta|>\epsilon)→0$となることである。
この条件を満足する場合、$\hat{\theta}_n$は一致推定量であるという。

基準3 漸近正規性

推定量が$X_1,\cdots,X_n$の複雑な関数であっても、中心極限定理により漸近分布は正規分布であることが多い。
漸近分布が正規分布である性質を漸近正規性と呼んだ。この性質を満たす推定量を漸近正規推定量と呼ぶ。
中心極限定理によれば、$\bar{X}$の漸近分布は、母集団分布に関係なく$N(\mu,\sigma^2\ /\ n)$であり、$\bar{X}$は漸近推定量である。

基準4 有効性

いかなる不偏推定量よりも分散が小さい推定量が存在すれば、それは非常に望ましい推定量であるといえる。
そのような推定量が存在するとき、有向推定量または最小分散不偏推定量と呼ぶ。
実際の推定では有向推定量を見つけることは難しいことが多い。
そこで、有効性の基準を少し緩めた漸近有効性と呼ばれる基準が使われることがある。
漸近有効性とは、漸近分布が正規分布おｔなる推定量のうち、その漸近分散が最小となる性質をいう。
$n→\infty$における有効性であり、この性質を満足する推定量を漸近的有向推定量と呼ぶ。

点推定の例

正規分布に関する推定

正規分布$N(\mu,\sigma^2)$の尤度関数は、

L(\mu,\sigma^2)=\prod_{i=1}^n\frac{1}{\sqrt{2\pi}\sigma}\exp\bigl(-(X_i-\mu)^2/2\sigma^2\bigr)\\
\log L(\mu,\sigma^2)=-n\log(\sqrt{2\pi}\sigma)-\sum_{i=1}^n(X_i-\mu)^2\ /\ 2\sigma^2

である。
$\mu,\sigma^2$で偏微分してそれぞれ0とおいて最尤推定量を求めると、

\bar{X}=\sum_{i=1}^nX_i\ /\ n,\hspace{5mm}S^2=\sum_{i=1}^n(X_i-\bar{X})^2\ /\ n

ただし、$S^2$は不偏でないので、$s^2=\sum_{i=1}^n(X_i-\bar{X})^2\ /\ (n-1)$を推定量とする方が一般的である。

二項分布に関する推定

母集団分布が標本$X_i=0,1$を出す母数$p$のベルヌーイ分布$Bi(1,p)$の二項分布であるとき、尤度関数の対数は、

$$
\log L(p)=\sum_{i=1}^n\bigl\{ X_i\log p + (1-X_i)\log (1-p) \bigr \}
$$

であるから、最尤推定量は$\hat{p}=\bar{X}$となる。

ポアソン分布に対する推定

母数$\lambda$の$Po(\lambda)$のとき、対数尤度関数は、

\log L(p)=\sum_{i=1}^n(-\lambda+X_i\log\lambda-\log X_i!)

である、したがって最尤推定量を求めると、

\hat{\lambda}=\bar{X}

となる。

一様分布に関する推定

母集団分布が区間$(a,b)$の一様分布($a<b$)であるとき、$a,b$が母数であるが、

\mu=(a+b)/2,\hspace{5mm}\sigma^2=(a-b)^2/12

である。これから$a,b$を解くと、

a=\mu-\sqrt{3}\sigma,\hspace{5mm}b=\mu+\sqrt{3}\sigma

であるから、モーメント法による推定量は、

a=\bar{X}-\sqrt{3}S,\hspace{5mm}b=\bar{X}+\sqrt{3}S

となる。
一方、最尤推定量は求め方は略すが、

$$
a=\min\{X_1,\cdots,X_n\}\hspace{5mm}b=\max\{X_1,\cdots,X_n \}
$$

区間推定

区間推定とは真の母数$\theta$が、ある区間$(L,U)$に入る確率を$1-\alpha$以上になるように保証する方法であり、

P\bigl(L\leq\theta\leq U\bigr)\geq 1-\alpha

となる確率変数$L,U$を求めるものである。
$L,U$はそれぞれ、下側信頼限界、上側信頼限界といわれる。
$1-\alpha$は信頼係数と呼び、区間$[L,U]$を$100(1-\alpha)%$信頼区間と呼ぶ。

正規母集団の母平均、母分散の区間推定

母平均の信頼区間

$\bar{X}$は正規分布$N(\mu,\sigma^2/n)$に従うから、これを標準化すると、

P\bigl(-Z_{\alpha/2}\leq \sqrt{n}(\bar{X}-\mu)\ /\ \sigma\leq Z_{\alpha/2}\bigr)=1-\alpha

これを$\mu$について解くと、

P\bigl(\bar{X}-Z_{\alpha/2}\sigma\ /\ \sqrt{n} \leq \mu \leq \bar{X}+Z_{\alpha/2}\sigma/\sqrt{n}\bigr)=1-\alpha

であり、$\mu$の信頼区間は

\bigl[\bar{X}-Z_{\alpha/2}\sigma\ /\ \sqrt{n},\bar{X}+Z_{\alpha/2}\sigma\ /\ \sqrt{n}\bigr]

となる。

また、$\sigma^2$が未知のときは標本分散$s^2$を使うと、$\sqrt{n}(\bar{X}-\mu)\ /\ s$は自由度$n-1$のt分布$t(n-1)$に従うから、

P\bigl(-t_{\alpha\ /\ 2}\leq \sqrt{n}(\bar{X}-\mu)\ /\ s\leq t_{\alpha\ /\ 2}\bigr)=1-\alpha

である。
これを$\mu$について解くと、

P\bigl(\bar{X}-t_{\alpha\ /\ 2}(n-1)s\ /\ \sqrt{n} \leq \mu \leq \bar{X}+t_{\alpha\ /\ 2}(n-1)s\ /\ \sqrt{n}\bigr)=1-\alpha

となる。したがって、母平均$\mu$の信頼区間$1-\alpha$の信頼区間は

\bigl[\bar{X}-t_{\alpha\ /\ 2}(n-1)s\ /\ \sqrt{n},\bar{X}+t_{\alpha\ /\ 2}(n-1)s\ /\ \sqrt{n}\bigr]

import pandas as pd
import numpy as np
import matplotlib.pyplot as plt
import japanize_matplotlib

d = pd.read_csv('pokemon_for_stats.csv').iloc[:,1:]
gen = ["第1世代"]*151+["第2世代"]*100+["第3世代"]*135+["第4世代"]*107+["第5世代"]*156+["第6世代"]*72+["第7世代"]*88+["第8世代"]*89
d['世代'] = gen

num_cols = ['重さ','高さ','HP','攻撃','防御','特攻','特防','素早さ','合計']
d.head()

ポケモンのデータの「攻撃」と「防御」を使って推定区間の計算を行う。

fig, ax = plt.subplots(1,2, figsize=(10,4))

ax[0].hist(d['攻撃'], density=True);
ax[0].set_title('攻撃')
ax[1].hist(d['防御'], density=True);
ax[1].set_title('防御')

Text(0.5, 1.0, '防御')

分散既知と分散未知の場合の区間推定の例を示す。
サンプル数を変化させ、どのように変化するか確認する。

from scipy.stats import norm, t

alpha = 0.1
# 分散既知の場合
sigma = np.std(d['攻撃'])

fig, ax = plt.subplots(1,3, figsize=(15,4))
for i, n in enumerate([4, 8, 40]):
    print('n=',str(n))
    # 各標本で平均と標準偏差の計算
    sample_mean = np.array([d['攻撃'].sample(n).mean() for _  in range(1000)])
    sample_std  = np.array([d['攻撃'].sample(n).std(ddof=1) for _  in range(1000)])

    # 推定値
    x_mean = np.mean(sample_mean)
    s_mean = np.mean(sample_std)

    # 正規分布のパーセント点(分散既知)
    z = norm.ppf(1-alpha/2)
    # t分布のパーセント点(分散未知)
    pt = t.ppf(1-alpha/2, n-1)

    print(' 信頼区間(分散既知): [',x_mean-z*sigma/np.sqrt(n),', ', x_mean+z*sigma/np.sqrt(n),']')
    print(' 信頼区間(分散未知): [',x_mean-pt*s_mean/np.sqrt(n),', ', x_mean+pt*s_mean/np.sqrt(n),']\n')

    # プロット
    ax[i].hist(sample_mean, color='lightblue');
    ax[i].vlines([x_mean-z*sigma/np.sqrt(n), x_mean+z*sigma/np.sqrt(n)], 0, 300, color='red', linestyles='--', label='分散既知');
    ax[i].vlines([x_mean-pt*s_mean/np.sqrt(n), x_mean+pt*s_mean/np.sqrt(n)], 0, 300, color='green', linestyles='--', label='分散未知');
    ax[i].set_xlim(20,140)
    ax[i].set_title('n='+str(n))
    ax[i].legend();

n= 4
 信頼区間(分散既知): [ 52.566587077420564 ,  101.33491292257943 ]
 信頼区間(分散未知): [ 44.1881177787593 ,  109.71338222124069 ]

n= 8
 信頼区間(分散既知): [ 59.15641804388649 ,  93.6408319561135 ]
 信頼区間(分散未知): [ 57.07162375784365 ,  95.72562624215634 ]

n= 40
 信頼区間(分散既知): [ 68.85280063280209 ,  84.27469936719791 ]
 信頼区間(分散未知): [ 68.66257308823583 ,  84.46492691176417 ]

分散既知の場合の方が、分散未知の場合より信頼区間が狭くなることが確認できる。
サンプル数の増加とともに区間は狭くなり、分散既知と分散未知の場合での差はなくなることがわかる。

母分散の信頼区間

$(n-1)s^2/\sigma^2$は自由度$n-1$の$\chi^2$分布$\chi^2(n-1)$に従うから、

P\bigl(\chi^2_{1-\alpha\ /\ 2}(n-1)\leq (n-1)s^2\ /\ \sigma^2\leq\chi^2_{\alpha\ /\ 2}(n-1)\bigr)=1-\alpha

となる。
これを$\sigma^2$について解くと、母分散$\sigma^2$の信頼係数$1-\alpha$の信頼区間は

\biggl[\frac{(n-1)s^2}{\chi^2_{\alpha\ /\ 2}(n-1)}, \frac{(n-1)s^2}{\chi^2_{1-\alpha\ /\ 2}(n-1)}\biggr]

となる。

ここでも、サンプル数を変化させながら標本分布及び信頼区間の変化の様子を確認する。

from scipy.stats import chi2
colors1 = ['lightblue', 'lightpink', 'lightgreen']
colors2 = ['blue', 'red', 'green']

for i, n in enumerate([4, 8, 40]):
    print('n=',str(n))
    # 標本の標準誤差の計算
    sample_std  = np.array([d['攻撃'].sample(n).std(ddof=1) for _  in range(1000)])

    # 推定値
    s_mean = np.mean(sample_std)

    # カイ二乗分布のパーセント点
    cl = chi2.ppf(1-alpha/2, n-1)
    cu = chi2.ppf(alpha/2, n-1)

    print(' 信頼区間: [',(n-1)*(s_mean**2)/cu,', ', (n-1)*(s_mean**2)/cl,']\n')

    # プロット
    plt.hist(sample_std**2, color=colors1[i], alpha=.5);
    plt.vlines([(n-1)*(s_mean**2)/cu, (n-1)*(s_mean**2)/cl], 0, 300, color=colors2[i], linestyles='--', label='n='+str(n));
plt.legend()

n= 4
 信頼区間: [ 6405.515842590535 ,  288.3986736841427 ]

n= 8
 信頼区間: [ 2614.165852083134 ,  402.7692865374226 ]

n= 40
 信頼区間: [ 1322.1435673159722 ,  622.5326786463075 ]

サンプル数の増加と共に分布の範囲は狭くなり、信頼区間も様くなることがわかる。

二つの正規母集団の母平均の差、母分散の比の区間推定

母平均の差の信頼区間

母集団分布$N(\mu_1,\sigma_1^2),N(\mu_2,\sigma_2^2)$である二つの正規母集団から個別に二つの標本$X_1,\cdots,X_m$と$Y_1,\cdots,Y_n$を抽出したときの、母平均の差$\mu_1-\mu_2$の区間推定について考える。
二つの母分散が等しく$\sigma_1^2=\sigma_2^2=\sigma^2$である場合、

s^2=\frac{1}{m+n-2}\biggl(\sum_{i=1}^m(X_i-\bar{X})^2+\sum_{j=1}^n(Y_j-\bar{Y})^2\biggr)

とすると、

t=\frac{(\bar{X}-\bar{Y})-(\mu_1-\mu_2)}{s\sqrt{\frac{1}{m}+\frac{1}{n}}}

は自由度$m+n-2$の$t$分布$t(m+n-2)$に従う。
したがって、

P\biggl(-t_{\alpha\ /\ 2}(m+n-2)\leq \frac{(\bar{X}-\bar{Y})-(\mu_1-\mu_2)}{s\sqrt{\frac{1}{m}+\frac{1}{n}}}\leq  t_{\alpha\ /\ 2}(m+n-2) \biggr)=1-\alpha

これを差$\mu_1-\mu_2$について解くと、信頼係数$1-\alpha$の信頼区間は

\biggl[\bar{X}-\bar{Y}-t_{\alpha\ /\ 2}(m+n-2)s\sqrt{\frac{1}{m}+\frac{1}{n}}, \bar{X}-\bar{Y}+t_{\alpha\ /\ 2}(m+n-2)s\sqrt{\frac{1}{m}+\frac{1}{n}} \biggr]

となる。
また、二つの母分散が等しいとと仮定できない場合は

s_1^2=\frac{1}{m-1}\sum(X_i-\bar{X})^2,\hspace{5mm}s_2^2=\frac{1}{n-1}\sum(Y_j-\bar{Y})^2

として、ウェルチの近似法を用いる、

t=\frac{(\bar{X}-\bar{Y})-(\mu_1-\mu_2)}{\sqrt{\frac{s_1^2}{m}+\frac{s_2^2}{n}}}

は自由度が

\nu=\frac{\bigl(\frac{s_1^2}{m}+\frac{s_2^2}{n} \bigr)^2}{\frac{(s_1^2/m)^2}{m-1}+\frac{(s_2^2/n)^2}{n-1}}

に最も近い整数$\nu^$の$t$分布$t(\nu^)$に近似的に従うから、信頼係数$1-\alpha$の信頼区間は、次式となる。

\biggl[\bar{X}-\bar{Y}-t_{\alpha\ /\ 2}(\nu^*)\sqrt{\frac{s_1^2}{m}+\frac{s_2^2}{n}}, \bar{X}-\bar{Y}+t_{\alpha\ /\ 2}(\nu^*)s\sqrt{\frac{s_1^2}{m}+\frac{s_2^2}{n}} \biggr]

例として、「攻撃」と「防御」の種族値(ゲーム内で個体差を生むための設定値)の平均の差の分布と信頼区間の確認をする。
それぞれのサンプル数は異なるものとし、数を増やして変化を見ることとする。
等分散であるとした場合と、等分散でない(非等分散)とした場合の両方の結果を確認する。

def calc_mean_var(df):
    return df.mean(),df.std(ddof=1)

# 標本のサイズ
M = [8,12,60]
N = [4,8,40]
alpha = 0.05

fig, ax = plt.subplots(1,2, figsize=(12,4))
for i, (m, n) in enumerate(zip(M, N)):
    print('n=',str(n),' m=',str(m))
    # 各標本の平均の計算
    sample_stats1 = np.array([calc_mean_var(d['攻撃'].sample(m)) for _  in range(1000)])
    sample_stats2 = np.array([calc_mean_var(d['防御'].sample(n)) for _  in range(1000)])
    
    x_bar = sample_stats1[:,0]
    y_bar = sample_stats2[:,0]
    std1 = sample_stats1[:,1]
    std2 = sample_stats2[:,1]
    
    # 各標本の平均値と標準偏差の推定値
    mu1 = np.mean(x_bar)
    mu2 = np.mean(y_bar)
    
    s1 = np.mean(std1)
    s2 = np.mean(std2)
    
    # 差の標準偏差
    s =  np.sqrt(((s1**2*(m-1))+(s2**2*(n-1)))/(m+n-2))
    
    # t統計量
    # 等分散
    t1 = ((x_bar - y_bar)-(mu1 - mu2))/(s*np.sqrt(1/m+1/n))
    # 非等分散
    t2 = ((x_bar - y_bar)-(mu1 - mu2))/(np.sqrt(s1**2/m+s2**2/n))

    # 非等分散の場合の自由度
    nu = int(((s1**2/m)+(s2**2/n))**2/((s1**2/m)**2/(m-1)+(s2**2/n)**2/(n-1)))
    
    # t分布
    x = np.arange(-3,3,0.02)
    t_dist1 = t.pdf(x, df=m+n-2)
    t_dist2 = t.pdf(x, df=nu)
    
    # パーセント点
    t10025 = t.ppf(1-alpha/2, m+n-2)
    t20025 = t.ppf(1-alpha/2, nu)
    
    print(' 信頼区間(等分散): [', mu1-mu2-t10025*s*np.sqrt(1/m+1/n),', ', mu1-mu2+t10025*s*np.sqrt(1/m+1/n),']')
    print(' 信頼区間(非等分散): [',mu1-mu2-t20025*np.sqrt(s1**2/m+s2**2/n),', ', mu1-mu2-t20025*np.sqrt(s1**2/m+s2**2/n),']\n')

    # プロット
    ax[0].hist((x_bar - y_bar), color=colors1[i], density=True, alpha=.5, bins=30, label='n='+str(n)+' m='+str(m));
    ax[1].hist(t1, color=colors1[i], density=True, alpha=.5, bins=30, label='n='+str(n)+' m='+str(m));
    ax[1].plot(x, t_dist1)
    ax[0].legend();
    ax[1].legend();

n= 4  m= 8
 信頼区間(等分散): [ -32.78010670567929 ,  44.19935670567929 ]
 信頼区間(非等分散): [ -35.432273021966175 ,  -35.432273021966175 ]

n= 8  m= 12
 信頼区間(等分散): [ -23.744287477621306 ,  31.8356208109546 ]
 信頼区間(非等分散): [ -23.968192767651328 ,  -23.968192767651328 ]

n= 40  m= 60
 信頼区間(等分散): [ -7.260140120298674 ,  16.47502345363202 ]
 信頼区間(非等分散): [ -7.236513827751887 ,  -7.236513827751887 ]

サンプル数が増えるとともに分布は狭くなり、信頼区間も狭くなることがわかる。
等分散であるとした場合の方が、非等分散であるとした場合より信頼区間は狭くなり、サンプル数が増加するとその差は小さくなることがわかる。
また、右図の分布を標準化したものは$t$分布に近いものとなった。

母分散の比の信頼区間

$(\sigma_2^2s_1^2)\ /\ (\sigma_1^2s_2^2)$は自由度$(m-1,n-1)$の$F$分布に従うから

P\bigl(F_{1-\alpha\ /\ 2}(m-1,n-1)\leq \frac{s_1^2}{s_2^2}\frac{\sigma_2^2}{\sigma_1^2} \leq F_{\alpha\ /\ 2}(m-1,n-1) \bigr)=1-\alpha

となる。
したがって、これを比$\sigma_2^2\ /\ \sigma_1^2$について解いて、この比の信頼係数$1-\alpha$の信頼区間を求めると

\bigl[F_{1-\alpha\ /\ 2}(m-1,n-1)s_2^2\ /\ s_1^2, F_{\alpha\ /\ 2}(m-1,n-1)s_2^2\ /\ s_1^2\bigr]

となる。

平均の差と同様に「攻撃」と「防御」の分散の比の信頼区間を調べる。

from scipy.stats import f

M = [8,12,60]
N = [4,8,40]
alpha = 0.05

for i, (m, n) in enumerate(zip(M, N)):
    print('n=',str(n))
    sample_stats1 = np.array([calc_mean_var(d['攻撃'].sample(m)) for _  in range(1000)])
    sample_stats2 = np.array([calc_mean_var(d['防御'].sample(n)) for _  in range(1000)])
  
    std1 = sample_stats1[:,1]
    std2 = sample_stats2[:,1]
    std1_std2 = std2**2 / std1**2
    
    s1 = np.mean(std1)
    s2 = np.mean(std2)
    s1_s2 =  s2**2 / s1**2
    
    F0025 = f.ppf(alpha/2, dfn=m-1, dfd=n-1)
    F0975 = f.ppf(1-alpha/2, dfn=m-1, dfd=n-1)
    print(' 信頼区間: [', F0025*s1_s2,', ', F0975*s1_s2,']')
    
    x = np.arange(0,15,0.2)
    f_dist = f.pdf(x, dfn=m-1, dfd=n-1)*s1_s2
    
    plt.hist(std1_std2, color=colors1[i], density=True, alpha=.5, bins=np.arange(0,15,0.1));
    plt.plot(x,f_dist, color=colors1[i], label='n='+str(n)+' m='+str(m))
    plt.xlim(0,6)
plt.legend()

n= 4
 信頼区間: [ 0.14543806848404814 ,  12.527314153379395 ]
n= 8
 信頼区間: [ 0.2550514643627952 ,  4.514714936050411 ]
n= 40
 信頼区間: [ 0.5599682893670989 ,  1.7835664159928286 ]

サンプル数が増加すると信頼区間が狭くなっていくことが確認できる。

二項、ポアソン母集団の各母数の信頼区間

母集団分布が母数$p$のベルヌーイ分布$Bi(1,p)$のとき、$\hat{p}=\bar{X}$で推定する。
二項分布$Bi(1,p)$に従う確率変数$X_i$は母平均$p$、母分散$p(1-p)$であるから、中心極限定理によって、$(\sum X_i-np)\ /\ \sqrt{np(1-p)}$の分布は$n$が大きい場合、標準正規分布で近似でき、

P\bigl(-Z_{\alpha\ /\ 2}\leq (\sum{X}_i-np)\ /\ \sqrt{np(1-p)} \leq Z_{\alpha\ /\ 2}\bigr)=1-\alpha

であり、$\hat{p}=\bar{X}=\sum X_i/n$として、$p$について解くと、

P\bigl(\hat{p}-Z_{\alpha\ /\ 2}\sqrt{p(1-p)\ /\ n}\leq p \leq \hat{p}+Z_{\alpha\ /\ 2}\sqrt{p(1-p)\ /\ n}\bigr)=1-\alpha

大数の法則によって$\hat{p}$は$p$の一致推定量であるから、$n$が大きい場合には、ほとんど$p$に等しいと考えられる。
したがって、

P\bigl(\hat{p}-Z_{\alpha\ /\ 2}\sqrt{\hat{p}(1-\hat{p})\ /\ n}\leq p \leq \hat{p}+Z_{\alpha\ /\ 2}\sqrt{\hat{p}(1-\hat{p})\ / \ n}\bigr)=1-\alpha

であり、有意水準$1-\alpha$の$p$の信頼区間は、近似的に

\bigl[\hat{p}-Z_{\alpha\ /\ 2}\sqrt{\hat{p}(1-\hat{p})\ /\ n}, \hat{p}+Z_{\alpha\ /\ 2}\sqrt{\hat{p}(1-\hat{p})\ /\ n}\bigr]

で求められる。

母集団分布が母数$\lambda$のポアソン分布$Po(\lambda)$である場合、平均、分散とも$\lambda$であるから、中心極限定理によって、$(\sum X_i-n\lambda)\ /\ \sqrt{n\lambda}$の分布は$n$が大きい場合、標準正規分布で近似できる。
したがって、$\lambda$の推定量$\hat{\lambda}=\bar{X}$を使って$\lambda$の信頼係数$1-\alpha$の信頼区間は

\biggl[\hat{\lambda}-Z_{\alpha\ /\ 2}\sqrt{\frac{\hat{\lambda}}{n}},\hat{\lambda}+Z_{\alpha\ /\ 2}\sqrt{\frac{\hat{\lambda}}{n}} \biggr]

で近似的に求めることができる。

ここでは、二項分布の平均値の信頼区間を確認する。
サンプル数を増やすことでどのように変化するかの確認をする。

from scipy.stats import binom
p = 0.2

fig, ax = plt.subplots(1,3, figsize=(12,4))
for i, n in enumerate([4,10,40]):
    print('n=',str(n))
    sample = binom.rvs(n=n,p=p,size=1000);
    bar_p = sample / n
    p_mean = sample.mean()/n
    
    x = np.arange(0,3,0.1)
    z = norm.pdf(x, loc=p_mean, scale=np.sqrt(p_mean*(1-p_mean)/n))
    b = binom.pmf(range(n), n=n, p=p)
    
    z_alpha = 1-norm.ppf(alpha/2)
    
    print(' 信頼区間: [', p_mean- z_alpha*np.sqrt(p_mean*(1-p_mean)/n),', ', p_mean+ z_alpha*np.sqrt(p_mean*(1-p_mean)/n),']')

    ax[i].hist(bar_p, bins=20, color=colors1[i], density=True, alpha=.5);
    ax[i].plot(x, z, color=colors1[i], label='norm')
    ax[i].plot(np.arange(n)/n, b*n, color=colors1[i], linestyle='--', label='binom')
    ax[i].set_xlim(0,1)
    ax[i].set_title('n='+str(n))
    ax[i].legend()

n= 4
 信頼区間: [ -0.39242968861788385 ,  0.8014296886178839 ]
n= 10
 信頼区間: [ -0.17511370028464124 ,  0.5703137002846412 ]
n= 40
 信頼区間: [ 0.01143660538238217 ,  0.38436339461761787 ]

サンプル数が増加すると平均値の分布は正規分布に近づき、近似的に信頼区間が求められることがわかる。

次回

仮説検定

参考書