６ 推定
最尤法とベイズ分析については、別のウェブサイトで扱っている。
ブートストラップ
本書の表6.5.1のデータに対して、p. 168のアルゴリズムによる相関係数のブートストラップによる分析をおこなうPythonサンプルスクリプトをリスト１に示す。スクリプトファイルは圧縮ファイルchap6_files.zipにまとめたので、ダウンロード・展開すれば自由に利用できる。
Jupyter Notebook用ファイルは、圧縮ファイルchap6_prg.zipとした。ダウンロード展開して自由に利用されたい。
リスト１のスクリプトを実行すると、ブートストラップによるサンプリングが始まり、終了すると相関係数rのサンプルのヒストグラムが図１のように表示される。
図１
図１において、モードはサンプルのKDE（kernel density estimation）により求められたものである。
図１のフォームを閉じると、端末には以下の出力が表示されている。
r = 0.47618832443367104
b = 1000 r = 0.5516597629675971
b = 2000 r = 0.3772372653251759
b = 3000 r = 0.549822390468271
b = 4000 r = 0.4537435776736825
b = 5000 r = 0.3652804340899794
b = 6000 r = 0.4982859805250741
b = 7000 r = 0.5062359139565838
b = 8000 r = 0.6079560522959311
b = 9000 r = 0.4939415586862759
b = 10000 r = 0.5437419717432709
mode = 0.4841100496363799
mean = 0.4746720291946079
median = 0.47734264439034907
95% interval = [0.310, 0.622]
最初のｒの値は、データに対して求められたそうかン係数である。それに続くｂとｒの値は、ブートストラップのサンプリングの途中における出力である。
最後の４行に、ブートストラップによるサンプリングのモード、平均値、中央値と95％区間が出力されている。
リスト１ ブートストラップのスクリプト例（s6_Bootstrap.py）
import numpy as np
import numpy.random as nr
import matplotlib.pyplot as plt
data = [
(44,
30), (44, 48), (60, 46), (61, 47), (40, 30),
(46,
44), (55, 66), (48, 51), (50, 55), (57, 62),
(42,
40), (59, 63), (51, 45), (53, 62), (54, 53),
(52,
30), (49, 46), (47, 54), (61, 46), (67, 67),
(53,
51), (46, 44), (41, 35), (58, 49), (67, 73),
(63,
61), (44, 63), (17, 42), (51, 54), (59, 52),
(45, 43), (60, 37), (51, 44), (55, 38), (49, 36),
(42,
42), (44, 48), (56, 52), (52, 47), (44, 65),
(61,
39), (47, 49), (30, 33), (43, 38), (48, 39),
(51,
47), (63, 61), (45, 57), (51, 44), (51, 51),
(43,
35), (44, 44), (49, 48), (36, 54), (60, 50),
(45,
44), (62, 62), (58, 62), (50, 63), (55, 58),
(56,
59), (63, 63), (32, 43), (48, 49), (39, 41),
(41,
54), (44, 50), (62, 71), (55, 63), (53, 54),
(42,
42), (55, 52), (66, 57), (36, 54), (50, 46),
(30,
46), (43, 49), (61, 48), (61, 53), (44, 45),
(44,
54), (45, 58), (54, 54), (57, 58), (52, 47),
(46,
43), (45, 50), (45, 38), (49, 67), (50, 38),
(34,
37), (47, 62), (29, 33), (54, 37), (53, 51),
(53,
50), (48, 51), (53, 49), (46, 43), (62, 54)
]
data = np.array(data)
N = len(data)
r = np.corrcoef(data, rowvar =
False)[0][1]
print('\nr =', r)
print()
#
#
ブートストラップ
#
B = 10_000
r_smpl = np.empty(B)
b = 1
while True:
idx = nr.randint(0, N, N)
smpl = data[idx]
r_smpl[b-1] =
np.corrcoef(smpl, rowvar = False)[0][1]
if b % 1000 == 0:
print('b =', b, ' r =',
r_smpl[b-1])
b += 1
if b > B:
break
def calc_mode(samples, a = 0.05,
n_points = 10_000):
"""
分布のサンプルからKDEによってモードを求める
"""
import numpy as np
import scipy.stats as ss
Lp, Up =
np.percentile(samples, [100 * a/2, 100 * (1 - a/2)])
coord = np.linspace(Lp, Up,
n_points)
est_pdf =
ss.gaussian_kde(samples).pdf(coord)
map_idx = np.argmax(est_pdf)
MAP_Est = coord[map_idx]
return MAP_Est, coord,
est_pdf, map_idx
#
# ブートストラップの結果の出力
#
r_mode = calc_mode(r_smpl) # モード
print('\nmode =', r_mode[0])
r_mean = np.mean(r_smpl) # 平均値
print('mean =', r_mean)
#
# 2.5, 50,
97.5パーセンタイル値
#
q025, q50, q975 =
np.percentile(r_smpl, [2.5, 50, 97.5])
print('median = ', q50)
print('95% interval = [{0:.3f},
{1:.3f}]'.format(q025, q975))
#
# ブートストラップのサンプルの」ヒストグラム
#
plt.hist(r_smpl)
plt.title('Mode(KDE) =
{0:.3f}'.format(r_mode[0]), fontsize = 16)
plt.xlabel('r', fontsize = 16)
plt.ylabel('Frequency', fontsize = 14)
plt.tight_layout()
plt.show()