統計学の備忘録的な学習メモ（1）に引き続き、Python（Scipyなど）を用いた統計学の勉強に関するメモです。（2）では主に様々な分布および検定についてまとめています。
個人的な備忘録の意味合いが強く、他記事と書き方が異なっており、説明なども少なめです。今後は予告なく内容を追加・修正するかもしれません。
※環境はJupyter Notebookの想定です。
インポートするライブラリ
import numpy as np
import pandas as pd
import scipy as sp
from scipy import stats
from scipy.special import gamma
from matplotlib import pyplot as plt
import seaborn as sns
sns.set()
二項分布のプロット
公正なコイン（きちんと表、裏がそれぞれ1/2の確率で出る）を50回投げて、そのうち表がk回出る確率をプロット。
n = 50
p = 1/2
k = np.arange(0, n + 1, 1)
pk = stats.binom.pmf(k, n, p)
plt.plot(k, pk, marker='o')
plt.xlabel("k")
plt.ylabel("P(X=k)")
plt.show()
二項分布を用いた検定
コインを20回投げたら表が15回出たとき、コインが公正と言えるかを検定する。
print(stats.binom_test(15, 20, 1/2)})
結果はおよそ0.0414となり、5%有意水準（0.05）では「コインが公正である」という帰無仮説は棄却される。
別の方法で確認する。表が5回以下となる確率、15 回以上となる確率をそれぞれ求める。当然ながら両者は一致し、またこれらの和が先程の0.0414と一致することが確認出来た。
lower = stats.binom.cdf(5, 20, 1/2)
upper = 1 - stats.binom.cdf(14, 20, 1/2)
print(f'lower = {lower:.4}')
print(f'upper = {upper:.4}')
print(f'upper + upper = {lower + upper:.4}')
lower = 0.02069
upper = 0.02069
upper + upper = 0.04139
t分布のプロット
さまざまな自由度（df）におけるt分布の形をプロットしてみる。標準正規分布も併せてプロットする。
df_array = np.arange(1, 10, 1)
fig = plt.figure(figsize=(10, 5))
ax = fig.add_subplot(111)
x = np.linspace(-4, 4, 100)
for i in range(0, len(df_array)):
y = stats.t.pdf(x, df=df_array[i])
ax.plot(x, y, label=f'df = {df_array[i]}', alpha=0.6)
ax.plot(x, stats.norm.pdf(x), label='norm', alpha=0.6)
ax.set_xlabel('t or z')
ax.set_ylabel('probability density')
plt.legend()
自由度の値が大きくなるほど、標準正規分布（norm表記）に近づくことが確認出来た。
t分布を用いた母平均の区間推定
あるサンプルについて、母集団の平均の95%信頼区間を推定する。
sample = np.array([10.1, 9.8, 11.1, 10.5, 9.6, 10.8])
mu = sp.mean(sample)
print(f'標本平均 = {mu:.4}')
sigma = sp.std(sample, ddof=1)
df = len(sample) - 1
print(f'自由度 = {df}')
se = sigma/sp.sqrt(len(sample))
interval = stats.t.interval(alpha = 0.95, df = df, loc = mu, scale = se)
print(f'下側信頼限界 = {interval[0]:.4}')
print(f'上側信頼限界 = {interval[1]:.4}')
標本平均 = 10.32
自由度 = 5
下側信頼限界 = 9.703
上側信頼限界 = 10.93
別の方法でも確認する。t.ppf()関数を用いて0.25%点、97.5%点を自分で計算し、信頼限界を求めてみる。
t_025 = stats.t.ppf(q = 0.025, df = df)
print(f'0.25%点 = {t_025:.4}')
t_975 = stats.t.ppf(q = 0.975, df = df)
print(f'97.5%点 = {t_975:.4}')
print(f'下側信頼限界 = {mu + t_025 * se:.4}')
print(f'上側信頼限界 = {mu + t_975 * se:.4}')
信頼区間は先程と確かに一致することが確認出来た。
0.25%点 = -2.571
97.5%点 = 2.571
下側信頼限界 = 9.703
上側信頼限界 = 10.93
※なお信頼区間の正しい解釈は中々難しい。今回の結果を「母平均の値は95%の確率で、9.703〜10.93の間に存在する」と考えるのは正しくない模様（母平均の真の値は知るのが困難だけで、あらかじめ確定した定数であり、信頼区間の間に存在する/しないは本来はっきり定まっているため、95%の確率で存在する…という考え方は出来ない）
正しい解釈は、また別のサンプルを母集団から同様の手順で抽出して、信頼区間の推定を行う…という手順を100回繰り返した場合、そのうち95回は信頼区間の間に母平均の値が存在する、という具合のようです。以下のサイトが参考となりました。
19-3. 95％信頼区間のもつ意味 | 統計学の時間 | 統計WEB
1群のt検定
例えば、平均厚み50 mm狙いで量産されている製品について、サンプルを10個抽出して厚みを実測したところ、以下の表のようになった（標本平均 50.67）。
| 50.2 |
51.4 |
49.7 |
51.1 |
49.4 |
50.9 |
52.1 |
51.4 |
50.9 |
49.6 |
この製品の平均厚みは本当に50 mmであると言えるか？。1群のt検定ではこうした問題に対する判断を行うことが出来る。有意水準は5%とする。
sample = np.array([50.2, 51.4, 49.7, 51.1, 49.4, 50.9, 52.1, 51.4, 50.9, 49.6])
mu = sp.mean(sample)
print(f'標本平均 = {mu:.4}')
df = len(sample) - 1
print(f'自由度 = {df}')
sigma = sp.std(sample, ddof=1)
se = sigma / sp.sqrt(len(sample))
t_value = (mu - 50)/se
print(f't値 = {t_value:.4}')
alpha = stats.t.cdf(t_value, df = df)
print(f'α = {alpha:.4}')
print(f'p値 = {(1 - alpha) * 2:.4}')
p値が約0.043となり、有意水準5%（0.05）を下回った。サンプルの平均厚みは50 mmとは有意に異なると言える。
標本平均 = 50.67
自由度 = 9
t値 = 2.349
α = 0.9783
p値 = 0.04337
1群のt検定は、専用の関数を用いてもっと簡単に行うことも可能である。
t_value, p_value = stats.ttest_1samp(sample, 50)
print(f't値 = {t_value:.4}')
print(f'p値 = {p_value:.4}')
t値、p値は確かに先程の例と一致している。
t値 = 2.349
p値 = 0.04337
母集団（すなわち製品全体）の平均厚みについて、信頼区間の推定も行ってみる。
interval = stats.t.interval(alpha = 0.95, df = df, loc = mu, scale = se)
print(f'下側信頼限界 = {interval[0]:.4}')
print(f'上側信頼限界 = {interval[1]:.4}')
信頼区間は50.02〜51.32となり、50を含まないことが確認出来た。
下側信頼限界 = 50.02
上側信頼限界 = 51.32
対応のある2群のt検定
例えば、5人の被験者（A〜E）について、異なる状態1、2における心拍数を測定した結果を測定した結果、以下の表のようになった。
| |
A |
B |
C |
D |
E |
| 状態1 |
105 |
99 |
94 |
88 |
110 |
| 状態2 |
96 |
93 |
87 |
90 |
102 |
状態1と2の違いが、心拍数に影響を及ぼしているだろうか？2群のt検定ではこうした問題に対する判断を行うことが出来る。有意水準は5%とする。
今回は同じ被験者に対するデータなので、状態1、2のデータはペアとして考えることが出来る。この場合は対応のある2群のt検定を行う。まず状態1、2の差分を取って、差分に対して平均が0となるかどうか、1群のt検定を行えば良い。
data = pd.DataFrame({
'person' : ['A','B','C','D','E','A','B','C','D','E'],
'state' : ['S1','S1','S1','S1','S1','S2','S2','S2','S2','S2'],
'heart_rate' : [105, 99, 94, 88, 110, 96, 93, 87, 90, 102],
})
state1 = np.array(data.query('state == "S1"')['heart_rate'])
state2 = np.array(data.query('state == "S2"')['heart_rate'])
diff = state1 - state2
t_value, p_value = stats.ttest_1samp(diff, 0)
print(f't値 = {t_value:.4}')
print(f'p値 = {p_value:.4}')
p値は約0.046となり、有意水準5%（0.05）を下回った。状態1と2の違いは、有意に心拍数に影響を及ぼしていると言える。
t値 = 2.85
p値 = 0.04638
対応のある2群のt検定も、専用の関数で簡単に行うことが可能である。結果は先程と全く同じとなる。
t_value, p_value = stats.ttest_rel(state1, state2)
print(f't値 = {t_value:.4}')
print(f'p値 = {p_value:.4}')
対応のない2群のt検定（Welchのt検定）
例えば、異なるクラスにおけるテストの平均点の比較、異なる工場で作られた同種製品の平均サイズの比較など、それぞれ独立したデータ同士には対応のない2群のt検定を適用する。
この場合、「かつては」2群の分散が等しいか異なるかをまず判定し（F検定で行う）、それによって2通りの方法を使い分ける必要がある、という考えが主流であったとのこと。しかし最近では…
- 分散が異なる場合用とされていた方法（Welchのt検定）は、分散が等しい場合に適用しても支障がないこと
- F検定 → t検定 の2段階を経ると、多重性の問題が生じること
といった理由によって、最初からWelchのt検定を行う考えが主流になっているとのことである。
【参考となった外部サイト】
等分散検定から t検定，ウェルチ検定，U検定への問題点
検定の多重性を分かりやすく解説します【F検定⇒t検定はダメ？】
24-4. 対応のない2標本t検定 | 統計学の時間 | 統計WEB
先程の心拍数のデータを、対応のないデータと仮定して、Welchのt検定を試してみる。
t_value, p_value = stats.ttest_ind(state1, state2, equal_var = False)
print(f't値 = {t_value:.4}')
print(f'p値 = {p_value:.4}')
使用データ自体は先程と同じだが、結果はもちろん異なっている。p値が約0.270となり有意水準5%（0.05）を上回っているので、有意差が得られないという結論となる。
t値 = 1.199
p値 = 0.2697
確認のために、Welchのt検定の式に従って自力でも計算してみる。state1、2の標本サイズを
、標本平均を
、不偏分散を
とすると、t値は以下の式で計算出来る。
自由度
は以下の式で近似的に求めることが出来る。
mean1 = state1.mean()
mean2 = state2.mean()
var1 = state1.var(ddof=1)
var2 = state2.var(ddof=1)
n1 = state1.size
n2 = state2.size
t_value = (mean1 - mean2) / np.sqrt((var1/n1 + var2/n2))
print(f't値= {t_value:.4}')
df = (var1/n1 + var2/n2)**2 / (var1**2/(n1**2*(n1 - 1)) + var2**2/(n2**2*(n2 - 1)))
alpha = stats.t.cdf(t_value, df = df)
print(f'p値 = {(1-alpha)*2:.4}')
結果は先程ときちんと一致することが確認出来た。
t値 = 1.199
p値 = 0.2697
さまざまな自由度（df）におけるカイ二乗分布の形をプロットしてみる。
df_array = np.arange(1, 10, 1)
fig = plt.figure(figsize=(10, 5))
ax = fig.add_subplot(111)
x = np.linspace(0, 10, 100)
for i in range(0, len(df_array)):
y = stats.chi2.pdf(x, df=df_array[i])
ax.plot(x, y, label=f'df = {df_array[i]}', alpha=0.6)
ax.set_xlabel('chi2')
ax.set_ylabel('probability density')
plt.legend()
例えば、ある製品を生産する工場AとBがあり、不良品の発生個数について以下のようなデータが得られたとする。
|
不良品 |
良品 |
| A |
780 |
37220 |
| B |
20 |
1980 |
工場Aの不良率は約2%、工場Bは約1%だが、これらは有意な差と言えるか？（適当なデータなので、どちらも不良品出し過ぎ！というツッコミは無しでw）
ピアソンのカイ二乗検定で検討する（有意水準は5%とする）。最初にデータのクロス集計表を作る。
data = pd.DataFrame({
'factory':['A', 'A', 'B', 'B'],
'quality':['good', 'defective', 'good', 'defective'],
'quantity':[37220, 780, 1980, 20]
})
cross = pd.pivot_table(
data=data,
values='quantity',
aggfunc=sum,
index='factory',
columns='quality'
)
cross
quality defective good
factory
A 780 37220
B 20 1980
集計表が出来れば、専用の関数で簡単に検定を行うことが可能である。
chi2, p_value, df, exp_freq_array = stats.chi2_contingency(cross, correction=False)
print(f'カイ二乗値 = {chi2}')
print(f'p値 = {p_value}')
print(f'自由度 = {df}')
print(f'期待度数表\n{exp_freq_array}')
p値は約0.001となり、5%（0.05）を下回るので、工場AとBの不良率には有意な差があると言える。
カイ二乗値 = 10.74
p値 = 0.001048
自由度 = 1
期待度数表
[[ 760. 37240.]
[ 40. 1960.]]
期待度数表について、手計算と比較して確認する。それぞれの合計値を記載した集計表を示す。
|
不良品 |
良品 |
合計 |
| A |
780 |
37220 |
38000 |
| B |
20 |
1980 |
2000 |
| 合計 |
800 |
39200 |
40000 |
もし工場A、Bの違いを無視してトータルで見た場合、不良率は800/40000 = 2%である。
このトータルの不良率通りになる場合、それぞれの工場における不良品個数は以下の通りになるはずである。これが期待度数となる。
- A工場：38000 * 2% = 760個
- B工場：2000 * 2% = 40個
良品も併せた期待度数表を以下に示す。確かに先程の出力と一致した。
|
不良品 |
良品 |
| A |
760 |
37240 |
| B |
40 |
1960 |
定義通りにピアソンのカイ二乗値も計算してみる。結果は10.74となり確かに先程の出力と一致した。
cross_array = np.array([[780, 37220],[20, 1980]])
exp_freq_array = np.array([[760, 37240],[40, 1960]])
chi2 = 0
for i in range(0, 2):
for j in range(0, 2):
chi2 += (cross_array[i][j] - exp_freq_array[i][j]) ** 2 / exp_freq_array[i][j]
print(f'カイ二乗値 = {chi2:.4}')
カイ二乗分布を用いて、母分散の区間推定を行ってみる。題材として東京大学教養学部統計学教室 編『基礎統計学I 統計学入門』東京大学出版会の練習問題11.6 ii)を解いてみる。東京における10日間の温度データから、その母分散
の信頼区間を信頼係数95%で求める。
data = np.array([21.8, 22.4, 22.7, 24.5, 25.9, 24.9, 24.8, 25.3, 25.2, 24.6])
df = data.size - 1
var = data.var(ddof=1)
chi2 = stats.chi2(df = df)
lower = df * var / chi2.ppf(1.0 - 0.05/2)
upper = df * var / chi2.ppf(0.05/2)
print(f'信頼区間 [{lower:.4}, {upper:.4}]')
得られた信頼区間は、書籍記載の練習問題の解答ときちんと一致した。
信頼区間 [0.9173, 6.462]
interval関数を用いた方法も試してみる。もちろん結果は先程と一致する。
data = np.array([21.8, 22.4, 22.7, 24.5, 25.9, 24.9, 24.8, 25.3, 25.2, 24.6])
df = data.size - 1
var = data.var(ddof=1)
interval = stats.chi2.interval(alpha = 0.95, df = df)
lower = df * var / interval[1]
upper = df * var / interval[0]
print(f'信頼区間 [{lower:.4}, {upper:.4}]')
等分散性の検定
正規分布に従う2つの群について、母集団の分散が等しいか否かを検証する場合はF検定を行う。ただしScipyのstatsモジュールには、F検定を簡単に直接行える関数は用意されていない模様なので、教科書などを参考に自前で実装してみる。
def ftest(d1, d2):
var1 = d1.var(ddof=1)
var2 = d2.var(ddof=1)
df1 = d1.size - 1
df2 = d2.size - 1
f_value = var1/var2 if var1 > var2 else var2/var1
p1 = stats.f.cdf(f_value, df1, df2)
p2 = stats.f.sf(f_value, df1, df2)
p_value = min(p1, p2) * 2
return f_value, p_value
np.random.seed(1)
data1 = stats.norm.rvs(loc = 50, scale = 5, size = 30)
data2 = stats.norm.rvs(loc = 50, scale = 10, size = 50)
f_value, p_value = ftest(data1, data2)
print(f" F-value = {f_value:.4}, P-value = {p_value:.4}")
結果は以下の通りであり、p値が有意水準5%（0.05）を下回ったので、実際に2群の分散が有意に異なることを検証できた。
F-value = 3.252, P-value = 0.0002677
参考までに、今回の2群をヒストグラムで描いて見ると下図のようになる。
なおstatsモジュールには、実はF検定とは別の、等分散性の検定を行う関数が用意されている（Bartlett検定、Levene検定）。これらは共に、本来は3群以上に対して適用する手法とのことであるが、2群でも適用可能であり、statsモジュールにおいてはF検定ではなくこれらを用いる想定である模様。
これらの原理や使い分けについてはあまり理解出来ていないが（F検定と同様、正規分布に従う2群であればBartlett検定を、そうでない場合はLevene検定を適用する？）、試しに先程F検定を行った2群に対して使用してみる。
data1 = stats.norm.rvs(loc = 50, scale = 5, size = 30)
data2 = stats.norm.rvs(loc = 50, scale = 10, size = 50)
f_value, p_value = ftest(data1, data2)
print(f"Ftest : F-value = {f_value:.4}, P-value = {p_value:.4}")
stat, p_value = stats.bartlett(data1, data2)
print(f"Bartlett : statistics = {stat:.4}, P-value = {p_value:.4}")
stat, p_value = stats.levene(data1, data2)
print(f"Levene : statistics = {stat:.4}, P-value = {p_value:.4}")
結果は以下の通り。いずれの手法でもp値は有意水準5%（0.05）を下回っている。
Ftest : F-value = 3.252, P-value = 0.0002677
Bartlett : statistics = 10.83, P-value = 0.0009991
Levene : statistics = 8.886, P-value = 0.003831
：浸漬後の重量
：溶媒Bの場合は1、そうでない場合は0となる変数
：溶媒Cの場合は1、そうでない場合は0となる変数
：誤差項
：パラメータ
