ブートストラッピング（平均値の分布）
サンプルに対する平均値の不安定さ（サンプル平均値の分布）をブートストラッピングによって推定する。すなわち、与えられたデータからランダムに繰り返しを許してサンプルサイズNのサンプリングを行い、そのサンプルに対して平均値を求める。このデータからのランダムサンプリングとそのサンプルに対する平均値の計算をn_simu回繰り返して平均値の分布を求める。Natsal-3のデータの平均値の分布のブートストラッピングによる推定を、Spiegelhalter(2019), p. xviに記載されているウェブサイトからダウンロードしたファイルを参考・使用したPythonスクリプトをリスト１のように作成した。入力ファイル02-4-sexual-partners-counts-x.csvは、Spiegelhalter(2019), p. xviに記載されているウェブサイトからダウンロードしたファイルに含まれているものである。
入力データファイルは、集計済みのものであるので、ブートストラッピングを適用するために、次のコードで個別データに展開している。
Data50 = []
TotalMen = 0
for P, C in zip(NumPartners50,
MenCount50):
if C > 0:
Data50 += [P] * C
TotalMen
+= C
Data50 = np.array(Data50)
# ndarray型に変換
個別データData50に対して、次のコードでブートストラッピングによる平均値の分布を求めている。
n_sim = 1000
# ブートストラッピングにおけるシミュレーション数
sample_sizes = [10, 50, 200, 760] # ブートストラッピングにおけるサンプルサイズ
for iboot in
range(4):
means = []
for isim
in range(n_sim):
sample = rng.choice(Data50, sample_sizes[iboot])
means.append(np.mean(sample))
リスト１のスクリプトを実行すると、図１のグラフが表示される。
図１
サンプルサイズN＝10、50、200、760に対する平均値の分布が描かれている。サンプルサイズが大きくなると、平均値の分布の幅は小さくなり、平均値が安定していくことが分かる。
import pandas as pd
import matplotlib.pyplot
as plt
import seaborn
as sb
import numpy
as np
import numpy.random
as nprd
# ダウンロードしたファイルの読み込み
d_csv = pd.read_csv('02-4-sexual-partners-counts-x.csv')
data = d_csv.values
NumPartners = d_csv['NumPartners']
MenCount = d_csv['MenCount']
#
# パートナー数が50以下の人に制限
#
NumPartners50 = NumPartners[NumPartners <= 50]
MenCount50 = MenCount[NumPartners <= 50]
#
# 集計データを個人データに展開
#
Data50 = []
TotalMen = 0
for P, C in zip(NumPartners50,
MenCount50):
if C > 0:
Data50 += [P] * C
TotalMen
+= C
Data50 = np.array(Data50)
# ndarray型に変換
rng = nprd.default_rng()
#
#
描画画面を４行１列の配置で４画面用意する
#
fig, ax = plt.subplots(nrows=4, ncols=1, figsize=(5,7), sharex=True, sharey=True)
n_sim = 1000
# ブートストラッピングにおけるシミュレーション数
sample_sizes = [10, 50, 200, 760] # ブートストラッピングにおけるサンプルサイズ
for iboot in
range(4):
means = []
for isim
in range(n_sim):
sample = rng.choice(Data50, sample_sizes[iboot])
means.append(np.mean(sample))
sb.kdeplot(means,
ax=ax[iboot])
ax[iboot].set_ylabel('Density')
ax[iboot].text(1,
1, f'N={sample_sizes[iboot]}', ha='left', va='top')
ax[3].set_xlabel('Number
of partners')
plt.tight_layout()
plt.show()
岡本安晴（2019） いまさら聞けないPythonでデータ分析．丸善出版
Spiegelhalter, D. (2019). The art of
statistics: Learning from data. Pelican