３ 検定
本書の第３章に関連するPythonスクリプト例として、２項分布、平均値の検定、ウェルチのｔ（Welch’s t）、等分散の検定、相関係数の検定、カイ２乗検定を行うものを以下に用意した。スクリプトファイルは圧縮ファイルchap3_files.zipにまとめた。ダウンロード・展開すれば自由に利用できる。
スクリプトファイルをJupyter Notebook用にまとめたものは、圧縮ファイルchap3_prgs.zipとして用意した。ダウンロード展開して自由に利用されたい。
なお、ｐ値についての議論があるが、これに関するウェブサイトも用意している。
ベイズ分析のPythonスクリプト例は、このウェブサイトなどに用意している。
２項分布（本書（3.1.1）式）による検定のスクリプトを、リスト3.1のように書いてみた。
リスト3.1 ２項分布による検定（s3_combi_test_code.py）
import scipy.special as sp
def cum_p_comb(p, N, k):
v = 0.0
for i in range(k+1):
v +=
sp.comb(N, i) * (p**i) * ((1-p)**(N-i))
return v
p = float(input('p = '))
N = int(input('N = '))
k = int(input('k = '))
print(cum_p_comb(p, N, k))
リスト3.1 を実行すると、例えば以下のようになる。
p = 0.5
N = 5
k = 3
0.8125
実行すると、確率ｐ、試行数Ｎ、成功試行数ｋの値の入力が求められるので、それぞれ設定していく。ｋの値を入力すると、計算が始まり、累積確率
が算出されて表示される。
リスト3.1のスクリプトをPythonの関数を利用するものに書き直すと、リスト3.1のようになる。
リスト3.２ ２項分布の関数を用いたスクリプト（s3_combi_test.py）
import scipy.stats as ss
p = float(input('p = '))
N = int(input('N = '))
k = int(input('k = '))
print(ss.binom_test(k, N, p,
alternative = 'less'))
実行は、リスト3.1の場合と同様である。
独立なデータの平均値の差の検定に用いられるｔ値を本書（3.2.3）式に従って求めるスクリプトをリスト3.3に示す。
リスト3.3 独立なデータの平均値の差のｔ値（s3_test_ind_code.py）
import numpy as np
import scipy.stats as ss
X = [5.128, 4.0, 4.587, 4.831]
Y = [4.082, 3.571, 3.891, 3.413]
print('X =\n', X)
print('Y =\n', Y)
mean_X = np.mean(X)
mean_Y = np.mean(Y)
M = len(X)
N = len(Y)
s_2 = (np.sum((X - mean_X)**2) +
np.sum((Y - mean_Y)**2)) / (M + N -2)
s_x_y_2 = (s_2/M) + (s_2/N)
t = (mean_X - mean_Y) / (s_x_y_2 **
0.5)
print('\nt =', t)
df = M + N -2
print('df =', df)
p = ss.t.cdf(t, df)
print('P(X < t) =', p)
print('P(X > t) =', 1-p)
リスト3.3を実行すると以下の出力を得る。
X =
[5.128, 4.0, 4.587, 4.831]
Y =
[4.082, 3.571, 3.891, 3.413]
t = 3.168566998316026
df = 6
P(X < t) = 0.9903230788463271
P(X > t) = 0.009676921153672913
ｔ値とともに、その左側の確率P(X < t)と右側の確率P(X > t)が出力されている。
リスト3.3のスクリプトをPythonの関数を用いたものに変更するとリスト3.4のようになる。
リスト3.4 独立なデータのｔ値をPythonの関数を用いて求めるスクリプト（s3_t_test_ind.py）
import numpy as np
import scipy.stats as ss
X = [5.128, 4.0, 4.587, 4.831]
Y = [4.082, 3.571, 3.891, 3.413]
print('X =\n', X)
print('Y =\n', Y)
M = len(X)
N = len(Y)
df = M + N -2
print('df =', df)
Results = ss.ttest_ind(X, Y,
alternative = 'less')
#print(Results)
#print('df =', df)
print('t =', Results[0])
print('P(X < t) =', Results[1])
print('P(X > t) =', 1 - Results[1])
リスト3.4を実行すると以下の出力を得る。
X =
[5.128, 4.0, 4.587, 4.831]
Y =
[4.082, 3.571, 3.891, 3.413]
df = 6
t = 3.168566998316026
P(X < t) = 0.9903230788463271
P(X > t) = 0.009676921153672913
ウェルチのｔを求めるプログラム例を拙著p.75の説明に基づき以下のように用意した。
import numpy as np
import scipy.stats as ss
X = [5.128, 4.0, 4.587, 4.831]
Y = [4.082, 3.571, 3.891, 3.413]
print('X =\n', X)
print('Y =\n', Y)
mean_X = np.mean(X)
mean_Y = np.mean(Y)
M = len(X)
N = len(Y)
sx_2 = np.sum((X - mean_X)**2) / (M -
1)
sy_2 = np.sum((Y - mean_Y)**2) / (N -
1)
s_x_y_2 = (sx_2/M) + (sy_2/N)
t = (mean_X - mean_Y) / (s_x_y_2 **
0.5)
print("\nt' =", t)
df = (((sx_2/M) + (sy_2/N))**2) / \
(((sx_2**2) / ((M**2) * (M-1))) + ((sy_2**2) / ((N**2) * (N-1))))
print("df' =", df)
p = ss.t.cdf(t, df)
print("P(X < t') =", p)
print("P(X > t') =", 1-p)
上記スクリプトを実行すると、以下の出力を得る。
X =
[5.128, 4.0, 4.587, 4.831]
Y =
[4.082, 3.571, 3.891, 3.413]
t' = 3.168566998316026
df' = 5.071464740943024
P(X < t') = 0.9878141600854486
P(X > t') = 0.012185839914551355
対応のあるデータの平均値の差の検定で用いるｔ値を、本書（3.2.6）式によって求めるスクリプトを、リスト3.5に示す。
リスト3.5 対応のあるデータの平均値の差の検定のｔ値を求めるスクリプト（s3_t_test_rel_code.py）
import numpy as np
import scipy.stats as ss
data = [[80, 85], [90, 88], [40, 60],
[55, 70], [30, 40]]
print(data)
d = np.array(data).T[0] -
np.array(data).T[1]
N = len(d)
d_mean = np.mean(d)
s_2 = np.sum((d - d_mean) ** 2) / (N -
1)
t = d_mean / ((s_2 / N) ** 0.5)
df = N - 1
print('\nt =', t)
print('df =', df)
p = ss.t.cdf(t, df)
print('P(X < t) =', p)
print('P(X > t) =', 1-p)
リスト3.5を実行すると、以下の出力が得られる。
X =
[80, 90, 40, 55, 30]
Y =
[85, 88, 60, 70, 40]
t = -2.5072881487467367
df = 4
P(X < t) = 0.03312441496721379
P(X > t) = 0.9668755850327863
リスト3.5をPythonの関数を用いたものに書き直すと、リスト3,6のようになる。
リスト3.6 Pythonの関数を用いた対応のあるデータのｔ値を求めるスクリプト（s3_t_test_rel.py）
import numpy as np
import scipy.stats as ss
data = [[80, 85], [90, 88], [40, 60],
[55, 70], [30, 40]]
print(data)
d = np.array(data).T[0] -
np.array(data).T[1]
N = len(d)
df = N - 1
Results = ss.ttest_1samp(d, popmean =
0, alternative = 'less')
print('\nt =', Results[0])
print('df =', df)
print('P(X < t) =', Results[1])
print('P(X > t) =', 1-Results[1])
リスト3.6を実行すると、以下の出力を得る。
X =
[80, 90, 40, 55, 30]
Y =
[85, 88, 60, 70, 40]
t = -2.5072881487467367
df = 4
P(X < t) = 0.03312441496721379
P(X > t) = 0.9668755850327863
等分散の検定を、本書（3.2.7）式により行うスクリプトを、リスト3.7に用意した。
リスト3.7 等分散の検定スクリプト（s3_F_test_EqVar.py）
import numpy as np
import scipy.stats as ss
X = [5.128, 4.0, 4.587, 4.831]
Y = [4.082, 3.571, 3.891, 3.413]
X = np.array(X)
Y = np.array(Y)
print('X =\n', X)
print('Y =\n', Y)
mean_X = np.mean(X)
mean_Y = np.mean(Y)
M = len(X)
N = len(Y)
s_X_2 = np.sum((X - mean_X)**2) / (M -
1)
s_Y_2 = np.sum((Y - mean_Y)**2) / (N -
1)
F = s_X_2 / s_Y_2
print('F =', F)
df1 = M - 1
df2 = N - 1
print('df1 =', df1)
print('df2 =', df2)
p = ss.f.cdf(F, df1, df2)
print('P(X < F) =',p)
print('P(X > F) =', 1-p)
リスト3.7を実行すると、以下の出力を得る。
X =
[5.128 4. 4.587 4.831]
Y =
[4.082 3.571 3.891 3.413]
F = 2.495833374961741
df1 = 3
df2 = 3
P(X < F) = 0.7638463585081164
P(X > F) = 0.23615364149188356
相関係数の検定を、本書（3.2.10）によって行うスクリプトを、リスト3.8に示す。
リスト3.8 相関係数の検定スクリプト（s3_corr_code.py）
import numpy as np
import scipy.stats as ss
data = [
(70,
168), (68, 165), (62, 159), (56, 156), (70, 170),
(64,
162), (59, 162), (70, 169), (68, 168), (69, 170),
(67,
168), (63, 166), (58, 163), (74, 170), (65, 168),
(60,
160), (65, 165), (72, 175), (61, 164), (77, 174),
(68,
164), (64, 161), (70, 166), (47, 147), (66, 165),
(68,
170), (73, 174), (74, 172), (63, 164), (83, 182),
(58,
159), (64, 166), (67, 167), (75, 174), (67, 163),
(68,
166), (68, 167), (66, 165), (63, 165), (55, 154),
(61,
163), (60, 160), (62, 164), (62, 160), (59, 160),
(51,
155), (59, 163), (68, 170), (56, 157), (79, 180)]
data = np.array(data)
print(data)
mx = np.mean(data.T[0])
my = np.mean(data.T[1])
SXY = np.sum((data.T[0] - mx) * (data.T[1] -
my))
SX2 = np.sum((data.T[0] - mx)**2)
SY2 = np.sum((data.T[1] - my)**2)
r = SXY / ((SX2 * SY2) ** 0.5)
print('\nr =', r)
N = len(data)
print('N =', N)
t = r * ((N - 2)**0.5) / ((1 -
r**2)**0.5)
print('t =', t)
df = N - 2
print('df =', df)
p = ss.t.cdf(t, df)
print('P(X < t) =', p)
print('P(X > t) =', 1-p)
リスト3.8を実行すると、以下の出力を得る。
r = 0.9428052005583538
N = 50
t = 19.595199469782504
df = 48
P(X < t) = 1.0
P(X > t) = 0.0
本書の表3.3.3のデータに対して式（3.3.1）によるカイ２乗検定を行うスクリプトを、リスト3.9のように作成した。
リスト3.9 適合度の検定スクリプト（s3_chi2fit_code.py）
import numpy as np
import scipy.stats as ss
pj = [1/2, 1/2]
Oj = [16, 4]
pj = np.array(pj)
Oj = np.array(Oj)
c = len(Oj)
Ej = np.empty(c)
print('pj =', pj)
print('Oj =', Oj)
N = np.sum(Oj)
print('N =', N)
for j in range(c):
Ej[j] = N * pj[j]
print('Ej =', Ej)
chi2 = 0
for j in range(c):
chi2 += ((Oj[j] - Ej[j])**2)
/ Ej[j]
print('ch2 =', chi2)
df = c - 1
print('df =', df)
p_chi2 = 1 - ss.chi2.cdf(chi2, df)
print('p =', p_chi2)
リスト3.9を実行すると、以下の出力を得る。
pj = [0.5 0.5]
Oj = [16 4]
N = 20
Ej = [10. 10.]
ch2 = 7.2
df = 1
p = 0.007290358091535665
リスト3.9のスクリプトをPythonの関数を用いたものに書き直すと、リスト3.10のようになる。
リスト3.10 Pythonの関数を用いた適合度検定スクリプト（s3_chi2fit.py）
import numpy as np
import scipy.stats as ss
pj = [1/2, 1/2]
Oj = [16, 4]
pj = np.array(pj)
Oj = np.array(Oj)
c = len(Oj)
Ej = np.empty(c)
print('pj =', pj)
print('Oj =', Oj)
N = np.sum(Oj)
print('N =', N)
for j in range(c):
Ej[j] = N * pj[j]
print('Ej =', Ej)
Results = ss.chisquare(Oj, Ej)
print('chi2 =', Results[0])
print('p =', Results[1])
リスト3.10を実行すると、以下の出力を得る。
pj = [0.5 0.5]
Oj = [16 4]
N = 20
Ej = [10. 10.]
chi2 = 7.2
p = 0.0072903580915356595
本書の表3.3.5のデータを、式（3.3.4）によって分布の違いを検定するスクリプトをリスト3.11のように作成した。
リスト3.11 分布の違いの検定スクリプト（s3_chi2_between_rows_code.py）
import numpy as np
import scipy.stats as ss
Oij = [[23, 27],
[38, 12]]
Oij = np.array(Oij)
print('Oij =\n', Oij)
r = len(Oij)
c = len(Oij[0])
print('r =', r)
print('c =', c)
ni = np.sum(Oij, axis=1)
print('ni =', ni)
Cj =np.sum(Oij, axis=0)
print('Cj =', Cj)
N = np.sum(Oij)
print('N =', N)
Eij = np.empty((r, c))
for i in range(r):
for j in range(c):
Eij[i][j] = ni[i] * Cj[j] / N
chi2 = 0.0
for i in range(r):
for j in range(c):
chi2
+= ((Oij[i][j] - Eij[i][j])**2) / Eij[i][j]
print('\nchi2 =', chi2)
df = (r - 1) * (c - 1)
p_value = 1.0 - ss.chi2.cdf(chi2, df)
print('p =', p_value)
print('df =', df)
リスト3.11を実行すると、以下の出力を得る。
Oij =
[[23 27]
[38 12]]
r = 2
c = 2
ni = [50 50]
Cj = [61 39]
N = 100
chi2 = 9.457755359394703
p = 0.002102582110668938
df = 1
本書の表3.3.6のデータを、式（3.3.7）によって独立性の検定を行うスクリプトをリスト3.12のように用意した。
リスト3.12 独立性の検定スクリプト（s3_chi2Ind_code.py）
import numpy as np
import scipy.stats as ss
Oij = [[29, 24],
[26, 21]]
Oij = np.array(Oij)
print('Oij =\n', Oij)
r = len(Oij)
c = len(Oij[0])
print('r =', r)
print('c =', c)
Ri = np.sum(Oij, axis=1)
print('ni =', Ri)
Cj =np.sum(Oij, axis=0)
print('Cj =', Cj)
N = np.sum(Oij)
print('N =', N)
Eij = np.empty((r, c))
for i in range(r):
for j in range(c):
Eij[i][j]
= Ri[i] * Cj[j] / N
print('Eij =\n', Eij)
chi2 = 0.0
for i in range(r):
for j in range(c):
chi2
+= ((Oij[i][j] - Eij[i][j])**2) / Eij[i][j]
print('\nchi2 =', chi2)
df = (r - 1) * (c - 1)
print('df =', df)
p_value = 1.0 - ss.chi2.cdf(chi2, df)
print('p =', p_value)
リスト3.12を実行すると以下の出力を得る。
Oij =
[[29 24]
[26 21]]
r = 2
c = 2
ni = [53 47]
Cj = [55 45]
N = 100
Eij =
[[29.15 23.85]
[25.85 21.15]]
chi2 = 0.0036495018429983615
df = 1
p = 0.951828212747989
リスト3.12をPythonの関数を用いたものに書き直すとリスト3.13のようになる。
リスト3.13 Pythonの関数を用いた独立性検定スクリプト（s3_chi2Ind.py）
import numpy as np
import scipy.stats as ss
Oij = [[29, 24],
[26, 21]]
Oij = np.array(Oij)
print('Oij =\n', Oij)
r = len(Oij)
c = len(Oij[0])
print('r =', r)
print('c =', c)
Results = ss.chi2_contingency(Oij,
correction = False)
print('Eij =\n', Results[3])
print('\nchi2 =', Results[0])
print('df =', Results[2])
print('p =', Results[1])
リスト3.13を実行すると以下の出力を得る。
Oij =
[[29 24]
[26 21]]
r = 2
c = 2
Eij =
[[29.15 23.85]
[25.85 21.15]]
chi2 = 0.0036495018429983615
df = 1
p = 0.951828212747989