２つの独立なサンプルのシミュレーション
データを母集団からのランダムなサンプリングと考えたとき、データはサンプリングごとにランダムに変動する。この変動の様子を見るためのシミュレーションプログラムを作成してみた。
２つの母集団からの独立なサンプリングについて考える。２つの独立なサンプルの平均値の検定について、岡本安晴「データを使いこなすための統計入門」で説明している。
平均がmu1、標準偏差がsgm1の正規分布からn1個のデータを採り、平均がmu2、標準偏差がsgm2の正規分からn2個のデータを採るとする。このときのデータの分布の様子を調べる。プログラム例をリスト１に示す。プログラムファイルは圧縮ファイルSimuTwoIndFiles.zipとした。ダウンロード展開すれば、自由に利用できる。
リスト１のプログラムを実行すると、母集団のパラメータとサンプルサイズの入力が求められるので、興味のある数値を設定していく。
mu_1 = 70
mu_2 = 80
sgm_1 = 25
sgm_2 = 20
n_1 = 50
n_2 = 60
サンプルサイズの設定まで終了すると、計算が始まる。
サンプリングは、次のスクリプトで行っている。
x1 = ss.norm.rvs(mu1, sgm1,
n1)
x2 = ss.norm.rvs(mu2, sgm2,
n2)
サンプルされたデータに対する統計量などは、次のスクリプトで求めている。
r = ss.ttest_ind(x1,
x2)
# Student
t = r.statistic
p = r.pvalue
vL, vU =
r.confidence_interval(0.95)
rw = ss.ttest_ind(x1, x2,
equal_var = False) # Welch
tw = rw.statistic
pw = rw.pvalue
vLw, vUw =
rw.confidence_interval(0.95)
サンプリングの結果などは、グラフに示される（図１）。
図１
タイトルの１行目と２行目に母集団のパラメータ値が表示されている。
３行目に、Studentのｔ値などが表示され、４行目にWelchのt’値などが表示されている。
グラフのウィンドウを閉じると、端末画面上に以下のように次のサンプリングシミュレーションに進むかどうか聞いてくる。
mu_1 = 70
mu_2 = 80
sgm_1 = 25
sgm_2 = 20
n_1 = 50
n_2 = 60
次に進む -> Nキー、 終了 -> Eキー
N or E ?
進む場合は、キーNを押す。以下に、Nキーを押してサンプリングを繰り返した例を示す。母集団のパラメータが同じでも、異なった結果が得られていることが分かる。
終了するときは、Eキーを押す。
リスト１ ２つの母集団からのサンプリング（simu_2popu.py）
import scipy.stats as ss
import numpy as np
import matplotlib.pyplot as plt
mu1 = float(input('mu_1 = '))
mu2 = float(input('mu_2 = '))
sgm1 = float(input('sgm_1 = '))
sgm2 = float(input('sgm_2 = '))
n1 = int(input('n_1 = '))
n2 = int(input('n_2 = '))
while True:
x1 = ss.norm.rvs(mu1, sgm1,
n1)
x2 = ss.norm.rvs(mu2, sgm2,
n2)
r = ss.ttest_ind(x1,
x2)
# Student
t = r.statistic
p = r.pvalue
vL, vU =
r.confidence_interval(0.95)
rw = ss.ttest_ind(x1, x2,
equal_var = False) # Welch
tw = rw.statistic
pw = rw.pvalue
vLw, vUw =
rw.confidence_interval(0.95)
str = f"mu1 = {mu1},
sgm1 = {sgm1}\nmu2 = {mu2}, sgm2 = {sgm2}\n" + \
f"Student: t = {t:.3f},
p = {p:.3f}" + \
f", 95%CI = [{vL:.1f},
{vU:.1f}]\n"+ \
f"Welch: t' =
{tw:.3f}, p' = {pw:.3f}" + \
f", 95%CI' =
[{vLw:.1f}, {vUw:.1f}]"
plt.hist(x1, alpha=0.5,
label='x1')
plt.hist(x2, alpha=0.5,
label='x2')
plt.title(str)
plt.legend()
plt.tight_layout()
plt.show()
print('次に進む -> Nキー、 終了 -> Eキー')
ck = input('N or E ? ')
if ck in ['E', 'e']:
break