区間推定1　クジ引き商品券の平均値

#点推定と区間推定
商品券のクジ引きを例に、点推定と区間推定を説明する。

ある商店街で商品券のクジ引きを行ったとしよう。景品には、600円券、300円券、100円券をそれぞれ大量に用意してあるが、クジ引き集計係の私はどの券がどの程度用意されているかを知らない。お客さんがクジ引きで当てた商品券の集計をした結果、くじ引き60回の平均金額は250円だったとする。この結果から、母集団の商品券の平均金額も同様に250円と推定するのが点推定である。点推定は理解しやすい推定だが、その推定値がどの程度信頼できるか分からない。そこで、点推定の値がどの程度の信頼度か推定する事を区間推定と言う。


#信頼区間(Confidence interval) 
上記の商品券のクジ引きで考えれば、未知の母数（母集団の特徴量）である商品券全ての金額平均を区間推定し、95％の確率で150円から250円の間になると推定できたとすると、信頼度95%（有意水準5%）の信頼区間は150～250円である、と表現する。
一般的な表現をすれば、未知の母数 θ が a ≦ θ ≦ b となる確率が(1-α)、すなわち　P[a≦ θ ≦b] = 1-α　となるように決定した [a, b] を、信頼度1-α(有意水準α)の信頼区間という。
統計の慣例として、信頼度は95%や99%（つまり有意水準は5%や1%）を用いる事が多い。
　

#平均の区間推定：商品券のクジ引きの場合
商店街のクジ引き商品券で、平均の区間推定を考える。

##平均の区間推定の背景・条件
くじ引き集計係である私は、60人の客のくじ引きを観察し、100円券が30枚、300円がで20枚、600円が10枚出た事を確認した。これら標本の平均金額は(100*30 + 300*20 + 600*10) / 60 = 250円である。商店街会長は、このクジ引きの商品券の平均金額は300円だと公表しているが、それは本当だろうか？私は母集団の平均値μ（クジ1枚で当たる商品券の平均金額）の信頼係数95%の信頼区間を推定して、会長の言葉の真偽を確かめようと考えた。


##平均の区間推定の計算
■標本数 $n$
引かれたクジの枚数nは次の通り。
n=10枚+20枚+30枚=60枚

■標本平均 $m$
引かれたクジの商品券金額は次の通り。
 $m=$(100*30 + 300*20 + 600*10) / 60 = 250円
大数の法則から、標本数を増やすと標本平均mは母集団平均μに近づくが、mとμの間には依然として差異がある状態である。
なお、標本平均は$\bar{x}$と表現する事も多い。

■標本分散$s^2$
標本分散は次の様になる。
　$標本分散s^2=\frac{1}{n}\sum_{i=1}^n (x_i-m)^2$

■[不偏分散][1] $u^2$　(unbiased variance)
標本数nが十分に大きくない場合、標本分散は母集団の分散より小さくなる。そこで期待値が母集団の分散に一致するように標本分散の算出式にn/(n-1)をかけたものが[不偏分散][1]である。
　$\begin{aligned}
　不偏分散u^2 &=\frac{1}{n-1}\sum_{i=1}^n (x_i-m)^2 \\\
    &=\frac{1}{59}((100-250)^2*30+(300-250)^2*20+(600-250)^2*10) \\\
    &=\frac{1}{59}(675000 + 50000 + 1225000)\\\
    &=33051
  \end{aligned}$

ところで、不偏分散は $u^2$や$v^2$と表現する事が多いが、標本分散で使われる$s^2$で表す事もある。nが大きければ不偏分散と標本分散はほぼ同じ値となるので、不偏分散の事を標本分散と表現して議論を進める事もある。nが100を超えるなら、標本分散だけ使用していても実用上問題になる事は少ないだろう。


■不偏標準偏差 $u$
不偏分散$u^2$の平方根から不偏標準偏差を求める。
　$u=\sqrt{33051} =181.8$

■標準誤差 $SE$
標準誤差（standard error）は推定量の標準偏差である。一般的に標準誤差とだけ言えば標本平均の標準偏差（standard error of mean）を意味する。
標準誤差の計算には中心極限定理が使われている。
中心極限定理は、平均μ、分散$σ^2$に従う母集団からサンプルサイズnの標本を抽出する時、nが大きくなるにつれて標本平均の分布は正規分布$N(μ,σ^2/n)$に近づくという定理である。つまり標本平均の分散は$σ^2/n$に近づき、標準偏差は$\sqrt{σ^2/n}$に近づくという定理である。

標準誤差SEは、標本平均の標準偏差であるから、中心極限定理より次の式が成り立つ。
　$SE=\sqrt{σ^2/n} =σ/\sqrt{n} $

本当に上記式が成り立っているだろうか。Pythonで確認してみる。

```Python:標準誤差
import numpy as np
import matplotlib.pyplot as plt

#標本平均を返す
def f_spl_mean(size,n_trial):
    spl_mean_arr = np.zeros(n_trial)
    for i in range(0,n_trial):
        sample = np.random.normal(loc=0, scale=2, size=size)
        spl_mean_arr[i] = np.mean(sample)
    return(spl_mean_arr)

#サンプルサイズの配列
size_arr = np.arange(start = 2, stop = 202, step = 2)
#標準誤差の配列
standard_error = 2 /size_arr**0.5
#標本平均の標準偏差の配列の準備
spl_m_std_arr  = np.zeros(len(size_arr))
np.random.seed(1)
for i in range(0,len(size_arr)):
    #母集団を平均0,標準偏差2の正規分布としてn個サンプリングし平均する処理を100回実施
    spl_mean = f_spl_mean(size_arr[i],100)
    #100回分の平均の不偏標準偏差
    spl_m_std_arr[i] = np.std(spl_mean,ddof = 1)

#グラフの描画
plt.figure(figsize=(6,4))
plt.title('Mean Standard')
plt.plot(size_arr,spl_m_std_arr,label='mean standard')
plt.plot(size_arr,standard_error,linestyle = 'dashed',linewidth=2,color='black',label='standard error')
plt.xlabel('sample size')
plt.ylabel('mean standard value')
plt.legend(loc='upper right',fontsize=12)
plt.grid(True)
plt.tight_layout()#一部が出力画像からはみ出る場合にグラフの位置サイズを調整
plt.savefig('figure.png')
plt.show()

```

・結果
![figure.png](https://qiita-image-store.s3.amazonaws.com/0/244452/c7294568-e57e-7516-59dc-f4a3e6087607.png)

上記プログラムでは、標準偏差が2（分散が4）の正規分布に従う母集団から標本をn個サンプリングして平均値を出す処理を100回繰り返し、100回の値から算出した不定標準偏差sを実線グラフで描画している。
また理論上の標本平均の標準偏差である標準誤差(2/√n)も破線グラフで描画している。
2つのグラフは一致しており、サンプルサイズnが大きくなるにしたがって平均の標準偏差(標準誤差)は小さくなっている。

今回のクジの問題では、私の立場では母分散$σ^2$が不明だが、標本から算出した不偏分散$u^2$が分かっている。不偏分散の期待値は母分散に一致するので、標準誤差SE（標本平均の標準偏差）は不偏標準偏差uを用いて次の様になる。
　$SE=u/\sqrt{n} = 181.8/7.746 = 23.47$


■標準正規分布に従うZ値
母分散の値が既知であればZ値を考えればよい。Z値は標準正規分布に従うので、そこから母平均μの区間推定が可能になる。
　$Z=(m-\mu)/SE=(m-\mu)/(σ/\sqrt{n})$

この式から、Z値とは標本平均mの標準化を表している事がわかるだろう。
正規分布を標準化して標準正規分布にする場合は、$(X-\mu)/\sigma\sim N(0,1)$という様に、確率変数Xから平均を引き標準偏差で割る。これと同様で、標本平均の標準化であるZ値は、標本平均から母平均を引き、標本平均の標準偏差である標準誤差で割るのである。

しかし現実の問題では、母分散が既知の場合は殆どない。母集団分散が既知の状況では、普通は母平均も既知なので、母平均を区間推定する事自体が無いのである。

今回のクジ引き商品券の例でも、当然母分散は不明である。よってZ値は計算できない。この様な場合はT値を考える。


■t分布に従うT値
T値はZ値と似ているが、不偏標準偏差で算出できる。
　$T=(m-\mu)/SE=(m-\mu)/(u/\sqrt{n})$

T値も標本平均mに対する標準化の計算を行っているが、おしい処で標準化に至っていない。たしかにT値の平均は0となるが、分散は1より大きくなる。これは、標本から計算された標準誤差を使用しているので、分散を膨らませてしまう為だ。つまりT値は標準正規分布には従わないという事である。
代わりに、T値は自由度n-1のt分布に従う。nは標本数を表すので、今回の例では自由度は59である。
自由度が100を超えるなら、t分布は標準正規分布にほぼ一致するので、実用上は近似的に標準正規分布表の値を利用しても良いだろう。しかし標本数を気にして使用する分布を切り替えるのも手間なので、どの様な標本数であろうとt分布を使用すれば間違いはない。
今回もt分布を使用して95%の信頼区間を求める。

■t分布を使用した区間推定
まずt分布表を参照し、自由度59・信頼度95%(有意水準0.05)のT値の臨界値を求めると2.001となるので、次の式が成り立つ。
※表計算ソフトの関数でもプログラムでもT値の臨界値は出力できるので、実務ではt分布表から臨界値を読み出す機会は少ないだろう。
　$P[-2.001<T<2.001]=0.95$

　$\begin{aligned}
　P[-2.001<T<2.001]
&=P[-2.001<(m-\mu)/(u/\sqrt{n})<2.001]\\\
&=P[m-2.001(u/\sqrt{n})<\mu< m+2.001(u/\sqrt{n})]\\\
&=P[250-2.001(181.8/\sqrt{60})<\mu< 250+2.001(181.8/\sqrt{60})]\\\
&=P[250-46.96<\mu< 250+46.96]\\\
&=P[203.04<\mu< 296.96]\\\
&=0.95
　\end{aligned}$


商品券の平均金額は、信頼度95%において203.04円から296.96円と推定する事ができた。会長が話していた平均300円は、推定した区間外にあるので、会長の嘘である可能性が高いと言えよう。


##T値による区間推定の前提条件
ｔ分布は正規分布する母集団から得た標本の平均に関する分布である。するとT値による区間推定を行う場合、母集団が正規分布である事が前提条件になるのだろうか。
結論から言えば、くじ引きの例の様に、ある程度標本数が大きければ、母集団が正規分布である必要は無い。正規分布ではない母集団から標本抽出した場合でも、ある程度標本数が大きければ、標本平均が中心極限定理により正規分布に従う。その為T値による区間推定が可能となるのである。

pythonで実験してみよう。
下記プログラムでは、次の2種の割合を出力して比較できる様にしている。
１．正規分布の標本から算出した信頼区間が母平均を含む割合
２．くじ賞金母集団の標本から算出した信頼区間が母平均を含む割合

くじ賞金の方は、母集団分布が正規分布とはかけ離れた形状ある事はご存知の通りである。この母集団分布でどの程度の精度で区間推定できるか、という趣旨で実験を行う。
プログラムの中ではくじ賞金の母集団や母平均を規定しているが、これは本来私が知り得ないデータで、区間推定の精度を確認する為に登場させたデータである。また、プログラムでは各割合を算出する為に10万回標本抽出して区間推定を計算しているが、これも区間推定の精度を確認する為のループ処理である。実務の区間推定では、1回だけ標本抽出をして、区間推定を行う事になる。

```python:区間推定の精度
import numpy as np
import matplotlib.pyplot as plt
from   scipy import stats

loop_count = 100000
n = 60 #サンプルサイズ
#ppf (Percent point function) t分布のt値を返す。
score = stats.t.ppf( 0.975, n -1) #片側0.975信頼区間(両側0.95と同様)、自由度n -1の値
print("T値{:.3f}".format(score))
np.random.seed(0)#同じ疑似乱数を生成する為、シード値を指定。

print("1.正規分布の標本から算出した信頼区間が母平均を含む割合")
t_high =np.array( [] )
t_low =np.array( [] )
for i in range(0,loop_count):
    #標準正規分布の母集団からのランダムなサンプル取得
    sample = np.random.normal(loc=0, scale=1, size=n)
    mu = np.mean(sample)
    std= np.std(sample, ddof = 1)
    se = std/ len(sample)**0.5
    t_high = np.append( t_high, mu + score * se )
    t_low  = np.append( t_low , mu - score * se )

#信頼区間が母平均を含む割合
print(1 - (len(t_high[t_high<0]) + len(t_low[t_low >0]) )/loop_count)


print("2.くじ賞金母集団の標本から算出した信頼区間が母平均を含む割合")
#くじ賞金母集団の配列ndarrayを作成　
arr_a=np.array([])
arr_a=np.hstack((arr_a,np.full(int(1),600)))#6の目が出たら600円
arr_a=np.hstack((arr_a,np.full(int(2),300)))#4,5の目が出たら300円
arr_a=np.hstack((arr_a,np.full(int(3),100)))#1,2,3の目が出たら100円
#母集団の平均
mu_population =np.mean(arr_a)

t_high =np.array( [] )
t_low =np.array( [] )
for i in range(0,loop_count):
    #母集団からのランダムにサンプル取得
    sample = np.random.choice(arr_a, n, replace=True)#replace=False抽出するデータに重複無し
    mu = np.mean(sample)
    std= np.std(sample, ddof = 1)
    se = std/ len(sample)**0.5
    t_high = np.append( t_high, mu + score * se )
    t_low  = np.append( t_low , mu - score * se )

#信頼区間が母平均を含む割合
print(1 - (len(t_high[t_high<mu_population]) + len(t_low[t_low >mu_population]) )/loop_count)
```
```:結果　サンプルサイズ60の場合
T値2.001
1.正規分布の標本から算出した信頼区間が母平均を含む割合
0.95021
2.くじ賞金母集団の標本から算出した信頼区間が母平均を含む割合
0.94769
```
1.正規分布の標本から算出した信頼区間が母平均を含む割合は、ほぼ95%となった。
2.くじ賞金母集団の標本から算出した信頼区間が母平均を含む割合は、母集団が左右非対称の分布なので、今回の標本数60程度だと多少95％を若干下回るが、実用上は問題ないレベルの精度と言えよう。

注意喚起の為に標本数が少ない場合も実験しておこう。プログラム上部のnを修正して実行する。

```:結果　サンプルサイズ3の場合
T値4.303
1.正規分布の標本から算出した信頼区間が母平均を含む割合
0.94961
2.くじ賞金母集団の標本から算出した信頼区間が母平均を含む割合
0.83464
```
1.正規分布の標本から算出した信頼区間が母平均を含む割合は、ほぼ95%となった。つまり、母数が正規分布なら標本数が小さくても推定可能である。推定可能ではあるが、標本数が小さければT値は大きくなり、実用的ではないので注意しよう。例えば「0円から500円の範囲と推定できた！」と主張しても、統計的に有意な結果なれど、実務的に有意義な結果ではない。
2.くじ賞金母集団の標本から算出した信頼区間が母平均を含む割合は、95%とかけ離れた値となった。標本数が小さすぎて、標本平均が正規分布にならない為である。

・結論
これら実験で、T値による区間推定において、標本数が大きければ母集団が正規分布でなくても良い事がわかるだろう。母集団が正規分布に従うなら、標本数が小さくてもT値による区間推定ができるとも言える。
なお、T値による区間推定だけでなく、Z値の区間推定でも同様の事が言える。



---
[区間推定2][2]

---

[1]:https://qiita.com/keymc021/items/8fcedc370f6563795573
[2]:https://qiita.com/keymc021/items/4f855790d9b2f7ee1bf6