目次
Python で統計学
母平均の差の検定
ここでは2つの正規母集団の差の検定を行ってみよう．具体的には100名のグループAと140名のグループBで得点の平均に差があるかどうかを検定したい状況を考えます．
データの確認
まずはライブラリを読み込んだあと，GitHub のリポジトリからグループ A とグループ B の2つのデータを読み込み，ヒストグラムを作成して眺めてみます．
ライブラリの読み込み
import numpy as np
import pandas as pd
import matplotlib.pyplot as plt
from scipy.stats import t
from scipy.stats import ttest_ind
from IPython.display import set_matplotlib_formats
# from matplotlib_inline.backend_inline import set_matplotlib_formats # バージョンによってはこちらを有効に
set_matplotlib_formats('retina')
グループ A のデータを読み込みます．
グループ A のデータの読み込み
urlA = "https://github.com/rinsaka/sample-data-sets/blob/master/groupA-scores.csv?raw=true"
dfA = pd.read_csv(urlA)
print(dfA)
ID score 0 101 237.1 1 102 207.2 2 103 207.9 3 104 202.8 4 105 196.7 .. ... ... 95 196 210.3 96 197 182.3 97 198 205.0 98 199 192.3 99 200 215.6 [100 rows x 2 columns]
グループ B のデータを読み込みます．
グループ B のデータの読み込み
urlB = "https://github.com/rinsaka/sample-data-sets/blob/master/groupB-scores.csv?raw=true"
dfB = pd.read_csv(urlB)
print(dfB)
ID score 0 501 190.5 1 502 159.8 2 503 160.2 3 504 208.4 4 505 227.4 .. ... ... 135 636 170.7 136 637 171.1 137 638 200.6 138 639 160.2 139 640 215.9 [140 rows x 2 columns]
グループ A の score を NumPy 配列 x に，グループ B の score を NumPy 配列 y にそれぞれ格納し，平均と分散を調べてみます．この結果，\(x\) の方が平均が大きく，分散が小さそうです．
x = dfA.loc[:, 'score' ].values
y = dfB.loc[:, 'score' ].values
print(f"xの平均： {np.mean(x):.3f}")
print(f"yの平均： {np.mean(y):.3f}")
print(f"xの分散： {np.var(x):.3f}")
print(f"yの分散： {np.var(y):.3f}")
xの平均： 201.403 yの平均： 196.688 xの分散： 276.624 yの分散： 517.243
2つの母集団でヒストグラムを作成してみよう．
x のヒストグラム
fig, ax = plt.subplots(1, 1, figsize=(6, 4))
ax.hist(x, bins=10)
ax.set_xlim(130, 270)
plt.show()
y のヒストグラム
fig, ax = plt.subplots(1, 1, figsize=(6, 4))
ax.hist(y, bins=10)
ax.set_xlim(130, 270)
plt.show()
2つのヒストグラムを眺めても，やはり \(x\) の平均がわずかに大きいように感じられる．また分散（ばらつき）は \(y\) の方が大きそうです．以降ではこの2つの正規母集団の母平均に差があるか（\(x\) の母平均の方が大きいか）について検定します．
検定の定式化
2つの母集団はそれぞれ正規分布 \(\mathcal{N}\left(\mu_1, \sigma_1^2\right)\), \(\mathcal{N}\left(\mu_2, \sigma_2^2\right)\) に従うものとします．なおここでは2つの母集団の母分散が未知である場合を考え，2つの母集団の母平均差があるかを検定します．帰無仮説は \(H_0~:~ \mu_1 - \mu_2 = 0\) を考え，対立仮説には右側検定の \(H_1~:~ \mu_1 - \mu_2 > 0\) を設定します．つまり，「母集団の母平均には差がない」という帰無仮説を棄却することで，\(x\) の方が平均が大きいことを示したい．
2つの母集団 \(x\), \(y\) のサンプル数をそれぞれ \(n_1\), \(n_2\) とします．このとき，ウェルチの近似法を用いると， \begin{eqnarray*} Z = \frac{\left(\bar{X} - \bar{Y}\right) - \left( \mu_1 - \mu_2 \right)} {\sqrt{U_1^2/n_1 + U_2^2/n_2}} \end{eqnarray*} は近似的に自由度が \begin{eqnarray*} \nu \approx \frac{\displaystyle \left(\frac{U_1^2}{n_1} + \frac{U_2^2}{n_2}\right)^2} {\displaystyle \frac{U_1^4}{n_1^2 \nu_1} + \frac{U_2^4}{n_2^2 \nu_2}} \end{eqnarray*} の \(t\) 分布に従うことになります．ここで， \(\bar{X}\) と \(\bar{Y}\) はそれぞれのサンプル平均であり，\(U_1^2\) と \(U_2^2\) は \begin{eqnarray*} U_1^2 &=& \frac{1}{n_1-1} \sum_{i=1}^{n_1} \left( X_i - \bar{X}\right)^2 \end{eqnarray*} \begin{eqnarray*} U_2^2 &=& \frac{1}{n_2-1} \sum_{i=1}^{n_2} \left( Y_i - \bar{Y}\right)^2 \end{eqnarray*} によって求めることができる不偏分散です．さらに，\(\nu_1 = n_1 -1\)，\(\nu_2 = n_2 -1\) です．また \(U_i^4 = (U_i^2)^2\) であることにも注意します． 対立仮説には \(H_1~:~ \mu_1 - \mu_2 > 0\) という右側検定を設定し，有意水準を \(\alpha\) とするとき，\(Z\) の実現値 \(z\) について，自由度 \(\nu\) の \(t\) 分布の上側 \(\alpha\) 点を \(p\) 値とします．このとき，
- \(p \leq \alpha\) なら，有意水準 \(\alpha\) で帰無仮説 \(H_0\) を棄却し，
- それ以外は，有意水準 \(\alpha\) で帰無仮説 \(H_0\) を棄却しない，
という結論が得られる．なお，対立仮説に \(H_1~:~ \mu_1 - \mu_2 \neq 0\) という両側検定を設定した場合には，\(t\) 分布の上側 \(\alpha/2\) 点を \(p\) 値とします．
Python で検定する
では Python で実際に検定を行ってみよう．上の \(Z\) の式に登場するデータ数，標本平均，不偏分散を順に求める．まずはデータ数 \(n_1\), \(n_2\) を求めて変数に格納します．
n_1 = x.size
n_2 = y.size
print(n_1, n_2)
100 140
次は標本平均を求める．
mean_1 = np.mean(x)
mean_2 = np.mean(y)
print(mean_1, mean_2)
201.40299999999996 196.68785714285715
さらに，不偏分散を求める．このとき ddof = 1 を指定する（つまり自由度を \(n-1\) とする）ことに注意します．詳細はここを参照してください．
u_1 = np.var(x, ddof=1)
u_2 = np.var(y, ddof=1)
print(u_1, u_2)
279.4178696969697 520.9642399794451
統計量 \(Z\) の実現値 \(z\) を計算します．
z = (mean_1 - mean_2)/np.sqrt(u_1/n_1 + u_2/n_2)
print(z)
1.8472510354620597
次は \(t\) 分布の自由度 \(\nu\) を計算します．なお，自由度は整数でなくても（実数であっても）定義できることに注意します．
nu = ((u_1 / n_1 + u_2 / n_2) ** 2) / ((u_1**2)/(n_1**2 * (n_1 - 1)) + (u_2**2)/(n_2**2 * (n_2 - 1)))
print(nu)
237.83722647173894
自由度 \(\nu\) の \(t\) 分布について，統計量 \(Z\) の実現値 \(z\) での \(\alpha\) を（つまり，\(p\) 値を）求める．
p = t.cdf(-np.abs(z), nu)
print('p_value : ', p)
p_value : 0.03297636855618921
したがって，有意水準を \(5\%\) とするとき，\(p = 0.03297636855618921\) が \(\alpha = 0.05\) よりも小さいことから「有意水準5%で帰無仮説 \(H_0\) を棄却する」という結論が得られた．
alpha = 0.05
if p <= alpha:
print("帰無仮説 H0 を棄却する")
else:
print("帰無仮説 H0 を棄却しない")
帰無仮説 H0 を棄却する
Scipy を使った検定
上で行ったような検定は Scipy ライブラリを用いてより簡単に実行することができます．つまり，scipy.stats の ttest_ind を使うだけで，上と同じ統計量 \(Z\) の実現値や \(p\) 値を得ることができます．具体的にはライブラリをインポートして，標本データを NumPy 配列に読み込んだあと次のコードを実行すると良いでしょう．
ttest_ind(x, y, equal_var = False, alternative = "greater")
Ttest_indResult(statistic=1.8472510354620597, pvalue=0.03297636855618921)
上の結果から \(p\) が pvalue=0.03297636855618921 として表示されています．これが \(5\%\) 以下であることから，帰無仮説を棄却するという結論を得ることができます．
なお，上の ttest_ind(x, y, equal_var = False, alternative = "greater") において，equal_var = False は2つの母集団の平均が異なることを指定しており，alternative = "greater" は \(H_1~:~ \mu_1 - \mu_2 > 0\) という右側検定を実行する指定を意味します．もしも \(H_1~:~ \mu_1 - \mu_2 < 0\) の左側検定ならば alternative = "less" を，\(H_1~:~ \mu_1 - \mu_2 \neq 0\) という両側検定ならば alternative = "two-sided" を指定すると良いでしょう．
また，ttest_ind( ) の結果はタプル形式で返されるので，次のように結果を変数に格納して利用しても良いでしょう．
statistic, pvalue = ttest_ind(x, y, equal_var = False, alternative = "greater")
print(statistic)
print(pvalue)
1.8472510354620597 0.03297636855618921
# 有意水準の設定
alpha = 0.05
if pvalue <= alpha:
print("帰無仮説 H0 を棄却する")
else:
print("帰無仮説 H0 を棄却しない")
帰無仮説 H0 を棄却する