７ 母集団を想定しない分析
本書で説明されているランダマイゼーション検定とサブサンプリングによる分析のPythonスクリプト例を以下のように用意した。なお、サンプルスクリプトのファイルは圧縮ファイルchap7_files.zipとしてまとめた。ダウンロード・展開すれば、自由に利用できる。
Jupyter Notebook用にまとめたものは圧縮ファイルchap7_prgs.zipとした。ダウンロード展開して自由に利用されたい。
本書の表7.1.1のデータに対してランダマイゼーション検定を行うPythonスクリプトをp. 173のアルゴリズムに従ってリスト１のように作成した。
リスト１を実行すると、以下の出力を得る。
[12 20 10 17 15 6 14 7
8 12]
diff_means = 5.4
I = 1000 N = 10000
I = 2000 N = 10000
I = 3000 N = 10000
I = 4000 N = 10000
I = 5000 N = 10000
I = 6000 N = 10000
I = 7000 N = 10000
I = 8000 N = 10000
I = 9000 N = 10000
I = 10000 N = 10000
337
p = 0.0337 se = 0.0018045583947326283
出力結果は、計算が乱数によるシミュレーションによっているので、実行ごとに確率的にランダムに変動するが、変動の大きさはseによって推測することができる。
リスト１ 本書の表7.1.1のランダマイゼーション検定のPythonスクリプト（s7_randomization.py）
import numpy as np
import numpy.random as nr
data = [12, 20, 10, 17, 15,
6,
14, 7, 8, 12]
data = np.array(data)
print(data)
diff_means = np.mean(data[:5]) -
np.mean(data[5:])
print('diff_means =', diff_means)
N = 10_000
# Random data
permutationの回数
smpls = np.empty(N)
smpls[0] = diff_means
Count = 1
I = 2
while True:
if I % 1000 == 0:
print('I =', I, ' N
=', N)
idx = nr.choice(10, 5,
replace = False)
s1 = []
s2 = []
for i in range(10):
if i
in idx:
s1.append(data[i])
else:
s2.append(data[i])
diff = np.mean(s1) -
np.mean(s2)
if diff >= diff_means:
Count += 1
smpls[I-1] = diff
I += 1
if I > N:
# Random data
permutationの回数がNに達したか？
break
print(Count)
p_hat = Count / N
se_hat = (p_hat * (1 - p_hat) / N) **
0.5
print('p =', p_hat, ' se =', se_hat)
本書の表7.2.3のデータのサブサンプルによる分析を行うPythonスクリプトを、p. 177のアルゴリズムに従って、リスト２のように作成した。
リスト２を実行すると、以下の出力を得る。
N = 50
N//2 = 25
(２組）ー（１組）:
平均 = 3.708 標準偏差 = 2.153
(３組）ー（１組）:
平均 = 3.720 標準偏差 = 3.725
(３組）ー（２組）:
平均 = 0.012 標準偏差 = 3.342
乱数によるシミュレーションであるので、出力にはサブサンプルによる分析を行うごとに確率的変動があることに注意。この確率的変動の大きさは、標準偏差の出力から推測できる。
リスト２ 本書の表7.2.3のサブサンプルによる分析のPythonスクリプト（s7_subsample.py）
import numpy as np
import numpy.random as nr
data1 = [
69, 55, 61,
53, 51, 50,
67, 47, 59,
69,
68, 55, 54,
73, 50, 51,
68, 62, 57,
64,
62, 74, 73,
83, 71, 35,
63, 70, 46,
73,
73, 64,
62, 77, 84,
54, 45, 63,
78, 52,
87, 68, 61,
56, 60, 53,
72, 42, 44,
53
]
data1 = np.array(data1)
data2 = [
55, 70, 55,
64, 77, 64,
61, 68, 66,
62,
74, 75, 66,
68, 56, 63,
75, 63, 59,
75,
72, 59, 65,
79, 49, 57,
72, 54, 61,
80,
62, 73, 55,
76, 77, 69,
64, 70, 53,
59,
64, 48, 65,
69, 68, 59,
65, 47, 60,
84
]
data2 = np.array(data2)
data3 = [
85, 84, 33,
94, 39, 32,
49, 40, 84,
47,
88, 49, 90,
83, 48, 44,
41, 42, 87,
45,
44, 49, 44,
46, 80, 85,
92, 48, 46,
83,
81, 87, 81,
83, 38, 80,
40, 90, 90,
92,
99, 91, 44,
89, 40, 47,
80, 49, 95,
34
]
data3 = np.array(data3)
N = len(data3)
print('N =', N)
print('N//2 =', N//2)
NSimu = 100
# シミュレーションの回数
Diff21 = np.empty(NSimu)
Diff31 = np.empty(NSimu)
Diff32 = np.empty(NSimu)
I = 1
while True:
subd1 = nr.choice(data1,
N//2, replace = False)
mean1 = np.mean(subd1)
subd2 = nr.choice(data2,
N//2, replace = False)
mean2 = np.mean(subd2)
subd3 = nr.choice(data3,
N//2, replace = False)
mean3 = np.mean(subd3)
Diff21[I-1] = mean2 - mean1
Diff31[I-1] = mean3 - mean1
Diff32[I-1] = mean3 - mean2
I += 1
if I > NSimu:
# シミュレーションの回数がNに達したか？
break
print()
print('(２組）ー（１組）:')
print('平均 =', f'{np.mean(Diff21):.3f}',
' 標準偏差 =', f'{np.std(Diff21):.3f}')
print('(３組）ー（１組）:')
print('平均 =', f'{np.mean(Diff31):.3f}',
' 標準偏差 =',
f'{np.std(Diff31):.3f}')
print('(３組）ー（２組）:')
print('平均 =', f'{np.mean(Diff32):.3f}',
' 標準偏差 =',
f'{np.std(Diff32):.3f}')