Pythonを使ってシミュレーションで学ぶ母平均の検定

# 目次
- 1.はじめに 
- 2.母平均の検定とは
    - 2.1 母平均の検定の目的
    - 2.2 t値の導入
    - 2.3 t値はt分布に従う
    - 2.4 判定ロジック
    - 2.5 有意水準に基づく棄却が出来ない場合何が判断出来るか

# 1.はじめに
　統計学を学びはじめて仮説検定につまづくというのは良くある話かとおもいます。かく言う私もそうです。
特に、確率変数を論理的に扱って帰無仮説を棄却する云々の仮説検定のロジックが腹落ちしませんでした。
本記事では、母平均の検定を題材に仮説検定の自分なりの理解を説明します。

Pythonの実行環境例：Anaconda3、Jupyter

# 2.母平均の検定とは
## 2.1 母平均の検定の目的
　母平均の検定は、品質管理などで良く使われます。例えば、1ロットで10000本の芯を製造する生産ラインがあったとします（あくまで例です）。ある検査工程では、同一ロットで製造した芯の平均（簡略化のために分散は考慮しません）が60mm以外ならば不合格とする場合、どのように判定するのが合理的でしょうか？全数調査では時間やコスト的に難しい場合、標本調査により合否を判定したいことがあると思います。
　ここで活躍するのが母平均の検定です。母平均では、標本を基に母平均が期待している値か否かを判定することを目的としています。

## 2.2 t値の導入
　母平均の検定にあたっては下記の検定統計量（以降、t値）を使用します。
$$t=\dfrac{(\hat{\mu}-\mu_0)}{\hat{\sigma}/N}$$
ここで$\hat{\mu}$は標本平均、$\hat{\sigma}$は標本分散、$N$は標本数であり、$\mu_0$が期待する母平均です。ではなぜ、わざわざ検定統計量を算出するのでしょうか？それは、「①標本が互いに独立で正規分布に従っている」と仮定する場合、t値が面白い性質を示すからです。以下、シミュレーションしながら性質を説明していきます。

## 2.3 t値はt分布に従う
　例えば、母集団が平均（以降、真の母平均）=60、分散=0.01の正規分布に従っていると仮定する場合、期待する母平均が60のt値の発生頻度は下記グラフのようになります。
![3.1.png](https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/35883/5726b78b-bb0a-73c3-abb2-e16db4d3d66d.png)
グラフ出力の実行コードは下記を参照。

```Python:t値の分布出力
import matplotlib
from scipy.stats import norm
import numpy as np
import matplotlib.pyplot as plt
#日本語表記に使用
matplotlib.rcParams['font.family']='Yu Mincho'

#母集団が正規分布に従うと仮定
#平均
mu=60
#分散
sigma=0.01
#標準偏差
stdev = np.sqrt(sigma)

#サンプル数
N=200

#期待する母平均
mu0=60

#試行回数
simu=10000
t_hist=[]

for i in range(simu):
    #標本平均・標本分散を算出する
    samples=norm.rvs(loc=mu, scale=stdev, size=N)
    # 標本平均⇒t検定の検定統計量算出に使用
    m=np.average(samples)
    # 標本標準偏差⇒t検定の検定統計量算出に使
    s=np.std(samples,ddof=1)
    #t検定の検定統計量
    t=(m-mu0)/(s/np.sqrt(N))
    t_hist=np.append(t_hist, t)

plt.hist(t_hist,bins=int(simu/100))
plt.title('母集団がN('+str(mu)+','+str(sigma)+')の正規分布に従う場合のt値の分布(期待する母平均='+str(mu0)+','+'標本数='+str(N)+')　試行回数='+str(simu)+'回')
plt.xlabel('t値')
plt.ylabel('発生頻度')
plt.xlim(-5,5)
plt.ylim(0,500)
plt.grid()
```

ご覧の通り正規分布のような形を示します。実はt値は、①の仮定の下では自由度がN-1のt分布（裾野の広い正規分布で自由度を上げると正規分布に近づく）に従います。**検定統計量が知られた確率分布に従うので、標本データから算出した検定統計量の発生確率を算出することが出来ます。**
![3.1-1.png](https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/35883/45845feb-0515-4758-4aea-8068929eded8.png)
　t値への理解を更に深めるために真の母平均を少し変化（$\pm 0.001$）させたときにt値の分布がどのように変わるかを見てみましょう。
![3.2.png](https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/35883/3fed38fc-5514-9c9f-138e-1068700a620d.png)
![3.3.png](https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/35883/89d2eaeb-80a9-5c89-f75e-1021e91a64aa.png)

このように、"真の母平均"="期待する母平均"の場合、t値の分布はt値=0を中心に左右対称になりますが、"真の母平均"<"期待する母平均"の場合は山のトップが負になり、"期待する母平均"<"真の母平均"の場合は山のトップが正になります。**真の母平均と期待する母平均の差により分布が左右にシフトします。**
　見やすくするためにt分布の出力結果を表示します。"真の母平均"="期待する母平均"の場合（青）と"期待する母平均"<"真の母平均"の場合（オレンジ）のt分布です。例えば、t値が2だったとします。この時、青の分布よりもオレンジの分布の方がt値が2になる確率が高いことが読み取れます。このことから直観的にですが、期待する母平均よりも真の母平均の方が大きい可能性が高いと解ります。
![3.4.png](https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/35883/9c56f8d6-9418-f473-02f8-b0192f66aecd.png)
グラフ出力の実行コードは下記を参照。

```Python:t分布出力
import numpy as np
import matplotlib.pyplot as plt
from scipy import stats
import matplotlib
matplotlib.rcParams['font.family']='Yu Mincho'

#自由度
Dof=200

#プロット数
plots=1000
tMatrix=np.zeros(NumberOfSamples)
tMatrix2=np.zeros(NumberOfSamples)
pMatrix=np.zeros(NumberOfSamples)

for i in range(plots):
    tMatrix[i]=stats.t.ppf(1-i/plots,Dof-1)
    tMatrix2[i]=talpMatrix[i]+1.5
    pMatrix[i]=stats.t.pdf(tMatrix[i],Dof-1)

plt.plot(tMatrix,pMatrix)
plt.plot(tMatrix2,pMatrix)
plt.title('自由度='+str(Dof)+'-1のt分布')
plt.ylabel('確率')
plt.xlabel('t値')
plt.ylim(0,0.5)
plt.xlim(-5,5)
plt.grid()
```
## 2.4 判定ロジック
　母平均の検定の目的は、"真の母平均"="期待する母平均"と判定できるかどうかでした。
例えば、200件の標本データを元に算出したt値が2の場合、どうすれば"真の母平均"="期待する母平均"がロジカルに判定できるでしょうか？論理学ではこのような場合、**背理法**を用いるかと思います。背理法では、ある仮定のもとで矛盾を導くことで仮定が成り立たないことを証明します。つまり、「"真の母平均"="期待する母平均"」という仮定のもとで矛盾を導くことが出来れば「"真の母平均"="期待する母平均"」が成り立たないことが言えます。仮説検定ではこの考え方を応用するために有意水準を導入します。**統計的手法を用いて判定したい仮定（帰無仮説という）が発生する確率が有意水準以下の場合、仮定が否定出来ること（統計の用語で棄却）を前提条件として設定**します。
　有意水準を導入すると判定ロジックは以下のようになります。

1. 有意水準を設定
2. 判定したい仮定（帰無仮説）を設定
3. 帰無仮説のもと検定統計量を算出
4. 3で算出した検定統計量から確率を算出
5. 4で算出した確率が1で設定した有意水準以下の場合、2で設定した帰無仮説を棄却⇒判定成功
6. 4で算出した確率が1で設定した有意水準より大きい場合、新たな帰無仮説を設定する

先述の例の場合、「"真の母平均"="期待する母平均"」という帰無仮説のもと$t\geq2$となる確率は約0.03（3%）です。有意水準として5%を設定した場合、帰無仮説が棄却できます。つまり、"真の母平均"$\neq$"期待する母平均"と判定できます。ここで、$t\geq2$としたのはt値の計算式よりt値が0から離れれば離れるほ真の母平均と期待する母平均が乖離する可能性が高いからです。

## 2.5 有意水準に基づく棄却が出来ない場合何が判断出来るか
　有意水準に基づく棄却が出来ない場合をもう少し説明します。これは、t分布のグラフをみただけではイメージが難しいです。t値が1の場合を考えます。**2.3**より、真の母平均と期待する母平均の差により分布が左右にシフトします。このとき次のグラフが役立ちます。真の母平均を変化させたとき$-1\leq t\leq 2$となる確率をシミュレーションしたグラフです。自由度が十分大きいt分布は$\pm2$に95.5%が分布するため、$t=1$をとりうる可能性が高い分布として$-1\leq t\leq 2$を指定しております。
![3.4-2-2.png](https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/35883/ece18f55-c0b6-e438-c78f-cec7fa9fc338.png)

$-1\leq t\leq 2$の場合、帰無仮説（グラフのx軸の0.000）における発生確率は約80%であるため有意水準5%と比較して大きく棄却できません。しかし、**帰無仮説が棄却できないからといって帰無仮説が正しいとは言えません。**グラフの通り、$-1\leq t\leq 2$となりうる真の母平均は多く存在します。従って、**帰無仮説が棄却出来ない場合は何も判断出来ません**

グラフ出力の実行コードは下記を参照。
※なお本コードは試行回数によっては処理に時間を要します。
そこで下記サイトを参考に進捗状況をコンソールに出力する処理を追加しております。


```Python:真の母平均を変化させたときの確率
import matplotlib
from scipy.stats import norm
import numpy as np
import matplotlib.pyplot as plt

#進捗表示用
import sys
import time
def progress(p, l):
    sys.stdout.write("\r%d / 100" %(int(p * 100 / (l - 1))))
    sys.stdout.flush()

#日本語表記に使用
matplotlib.rcParams['font.family']='Yu Mincho'

#母集団が正規分布に従うと仮定
_mu=60
sigma=0.01
stdev = np.sqrt(sigma)

#サンプル数
N=200

#期待する母平均
mu0=60

#試行回数
simu=1000
t_hist=[]

# 母平均の変更
mu_sim =100
delta=stdev/mu_sim
mu_start=_mu-stdev/2
mu_hist =np.zeros(mu_sim)
t_hist_geq2=np.zeros(mu_sim)
_t=-1
_t2=2

cntAll=0
for j in range(mu_sim):
    progress(j, mu_sim)
    mu=mu_start+j*delta
    mu_hist[j]=mu-mu0
    cnt=0

    for i in range(simu):
        #標本平均・標本分散を算出する
        samples=norm.rvs(loc=mu, scale=stdev, size=N)

        # 標本平均⇒t検定の検定統計量算出に使用
        m=np.average(samples)

        # 標本標準偏差⇒t検定の検定統計量算出に使
        s=np.std(samples,ddof=1)

        #t検定の検定統計量
        t=(m-mu0)/(s/np.sqrt(N))
        if(t>=_t and t<=_t2):
            cnt=cnt+1

    t_hist_geq2[j]=cnt

plt.plot(mu_hist,t_hist_geq2/simu)
plt.title('真の母平均を変化させたときの'+str(_t)+'<=t<='+str(_t2)+'となる確率（期待する母平均='+str(mu0)+',真の母分散='+str(sigma)+')　試行回数='+str(simu)+'回')
plt.xlabel('真の母平均-期待する母平均')
plt.ylabel(str(_t)+'<=t<='+str(_t2)+'となる確率')
plt.xlim(-stdev/2,stdev/2)
plt.grid()
```