独立なサンプルの平均値の差のベイズ分析
２つの独立なサンプルの平均値の差についてベイズ分析を行う。
A市から１５人、B市から２０人をランダムに選び血圧を測ったところ図１に示す結果を得たとする。
図１ data.xlsx
データの形式は、１行目に変数名を書き、２行目以降にデータ値が書かれている。第１列は各データの識別用で数値でも文字列でもよい。２列目に第１グループ（図１の場合はA市）のデータを書き、３列目に第２グループ（図１の場合はB市）のデータを書く。人数が揃わないときは、少ない方のデータは欠損値として空欄にしておく。
図１の形式のデータに対して以下の確率モデルを設定して平均値の差のベイズ分析を行う。
上のモデルに対するStanスクリプト例をリスト１としてこのウェブサイトの後半に示す。なお、サンプルファイルは圧縮ファイルfilesBysDifMuInd.zipとしてまとめた。ダウンロード解凍すれば自由に利用できる。
リスト１のStanスクリプトを利用するPythonスクリプトをリスト２としてこのウェブサイトの後半に示す。
リスト２のスクリプトを実行すると、入力データファイル名の設定を求めるダイアログウィンドウが表示される（図２）。
図２
入力データファイルを設定する。図２では、図１のファイル名data.xlsxが設定されている。データファイル名を設定して「開く」ボタンをクリックすると、計算結果などが書き出される出力ファイル名の設定を求めるダイアログウィンドウが表示される（図３）。
図３
適当なテキストファイル名を設定する。ファイル名を設定して「保存」ボタンをクリックすると、グラフ画像を保存するファイルの名前の設定を求めるダイアログウィンドウが表示される（図４）。
図４
画像は複数枚保存されるので、図４のダイアログウィンドウで設定した名前を基に複数の名前が作成されて保存される。設定した名前が図４のように「MyFig.png」であれば、「MyFig_XXX.png」という形式の名前（XXXの部分は図によって変わる）で画像が保存される。保存ファイル名を設定して「保存」ボタンをクリックすると、Stanスクリプトのコンパイルが始まり、コードのビルド後、MCMCによるサンプリングが始まる。サンプリングが終了すると、まず平均値パラメータの事後分布が表示される（図５）。
図５
図５のウィンドウを閉じると、平均値の差の事後分布が表示される（図６）。
図６
差の事後分布の95％CIが[-11.77, -6.11]と負の領域にあるので、A市とB市の平均値の差は負である、すなわちA市の平均値の方が低いことが分かる。
図７のウィンドウを閉じると、標準偏差の比の事後分布が表示される（図７）。
図７
標準偏差の比の事後分布が、平均値が0.957と１に近く、95％CIも[0.556, 1.584]と１を内部の中央付近に含む位置である。標準偏差は等しいと考えてよいようである。
図７のウィンドウを閉じると、２つの標準偏差それぞれの事後分布が表示される（図８）。２つの事後分布の重なりが明瞭である。
図８
図８のウィンドウを閉じると、図９の事後予測分布のグラフが表示される。
図９
２つの群（A市とB市）からそれぞれランダムにサンプルを採り、その差を求めたものである。２つの群は、パラメータ値が事後分布から取られた値である正規分布に従うものとしてシムレーションによって求めている。A市とB市からランダムに１人ずつ選んで比べたとき、A市の人の方が値が高いのは100回中数回のレベルであるという予測である。
図９のウィンドウを閉じると、プログラムの実行終了である。実行終了後、図３で指定した出力ファイルを開くと以下のようになっている。PythonのIDLEでは開けないことがあるが、そのときはメモ帳やワードパッドなど他のエディタで開けばよい。
入力データファイル:
D:/yasuharu/www/LaCoocan/Python/IntroPyStat/BysDiffMeanInd/filesBysDifMuInd/data.xlsx
A市
[126, 122, 120, 132, 129, 124, 124,
125, 124, 124, 121, 121, 130, 125, 129]
B市
[126, 139, 137, 130, 139, 132, 128,
137, 134, 137, 131, 133, 136, 131, 130, 134, 141, 135, 136, 133]
Inference for Stan model:
anon_model_ccbd0563488bdd44c9dc307138ff4420.
4 chains, each with iter=2000;
warmup=1000; thin=1;
post-warmup draws per chain=1000,
total post-warmup draws=4000.
mean
se_mean sd 2.5% 25% 50% 75% 97.5% n_eff Rhat
muA 125.05 0.02 1.06 122.99 124.38 125.02 125.72
127.19 3238 1.0
muB 133.97 0.02 0.97 132.13 133.32 133.96 134.6 135.94 3642 1.0
sgmA 3.91 0.02 0.83 2.68 3.33 3.78 4.34 5.94 2665 1.0
sgmB 4.21 0.01 0.75 3.03 3.66 4.12 4.64 5.96 2913 1.0
mu_d -8.92 0.02 1.43 -11.77 -9.86 -8.9 -7.98 -6.1 3292 1.0
sgm_r 0.96 4.8e-3 0.26 0.56 0.77 0.92 1.1 1.59 2975 1.0
lp__ -62.36 0.04 1.5 -66.05 -63.09 -62.02
-61.24 -60.48 1639 1.0
Samples were drawn using NUTS at Wed
Jul 8 10:21:37 2020.
For each parameter, n_eff is a crude
measure of effective sample size,
and Rhat is the potential scale
reduction factor on split chains (at
convergence, Rhat=1).
muA:
mean = 125.052
95%CI = [122.998, 127.192]
Q1 = 124.377, Median = 125.022, Q3 = 125.719
muB:
mean = 133.969
95%CI = [132.134, 135.935]
Q1 = 133.315, Median = 133.956, Q3 = 134.595
mu_d:
mean = -8.917
95%CI = [-11.772, -6.111]
Q1 = -9.856, Median = -8.899, Q3 = -7.984
sgm_r:
mean = 0.957
95%CI = [0.556, 1.584]
Q1 = 0.770, Median = 0.917, Q3 = 1.102
sgmA:
mean = 3.911
95%CI = [2.682, 5.938]
Q1 = 3.331, Median = 3.780, Q3 = 4.344
sgmB:
mean = 4.214
95%CI = [3.033, 5.954]
Q1 = 3.663, Median = 4.120, Q3 = 4.640
リスト１ Stanスクリプト（IndDiffMu2Sgm.stan）
data {
int NA;
int NB;
real DataA[NA];
real DataB[NB];
}
parameters {
real muA;
real muB;
real<lower = 0.0>
sgmA;
real<lower = 0.0>
sgmB;
}
transformed parameters {
real mu_d;
real sgm_r;
mu_d = muA - muB;
sgm_r = sgmA / sgmB;
}
model {
muA ~ normal(0.0, 1000.0);
muB ~ normal(0.0, 1000.0);
sgmA ~ exponential(0.001);
sgmB ~ exponential(0.001);
for (i in 1:NA)
DataA[i] ~ normal(muA, sgmA);
for (i in 1:NB)
DataB[i] ~ normal(muB, sgmB);
}
リスト２ Stanスクリプトを用いるPythonスクリプト（AnalIndMeans2SgmBys.py）
import tkinter as tk
import tkinter.filedialog as fd
import pandas as pd
import numpy as np
import pystan
import matplotlib.pyplot as plt
import seaborn as sb
import scipy.stats as ss
def askinfilename(msg = None, types =
[('', '*.*')]):
"""
入力用ファイル名の設定
"""
rt = tk.Tk()
rt.withdraw()
flnm = fd.askopenfilename(title
= msg, filetypes = types)
rt.destroy()
return flnm
def askoutfilename(msg = None, types =
[('', '*.*')], deftype = '.txt'):
"""
出力用ファイル名の設定
"""
rt = tk.Tk()
rt.withdraw()
flnm = fd.asksaveasfilename(title
= msg, filetypes = types,
defaultextension = deftype)
rt.destroy()
return flnm
finnm = askinfilename('入力データExcelファイル', types = [('Excel(*.xlsx)', '*.xlsx')])
txtoutnm = askoutfilename(msg = '出力ファイル', types = [('Text
(*.txt)', '*.txt')])
txtout = open(txtoutnm, 'w')
foutnm = askoutfilename(msg = '画像保存ファイル', types = [('画像 (*.png)', '*.png')],
deftype = '.png')
xlsx_fl = pd.ExcelFile(finnm)
data_xlsx =
pd.read_excel(xlsx_fl, header =
None)
data = data_xlsx.values
print(data)
txtout.write('\n入力データファイル:\n')
txtout.write(' {0}\n'.format(finnm))
VarNameA = data[0][1]
VarNameB = data[0][2]
dataA = []
dataB = []
for i in range(1,len(data)):
if not
pd.isnull(data[i][1]):
dataA.append(data[i][1])
if not
pd.isnull(data[i][2]):
dataB.append(data[i][2])
txtout.write('\n{0}\n'.format(VarNameA))
txtout.write('{0}\n'.format(dataA))
txtout.write('\n{0}\n'.format(VarNameB))
txtout.write('{0}\n'.format(dataB))
StanData = {'NA': len(dataA), 'DataA':
dataA,
'NB': len(dataB), 'DataB': dataB}
sm = pystan.StanModel(file =
'IndDiffMu2Sgm.stan')
fit = sm.sampling(data = StanData,
n_jobs = 1)
print(fit)
txtout.write('\n{0}\n'.format(fit))
# 日本語表記の設定
plt.rcParams['font.family'] = ['Yu
Gothic', 'Noto Sans CJK JP', 'Hiragino Maru Gothic Pro']
muA_p025, muA_p25, muA_p50, muA_p75,
muA_p975 = \
np.percentile(fit['muA'],
[2.5, 25, 50, 75, 97.5])
muA_mean = fit['muA'].mean()
muB_p025, muB_p25, muB_p50, muB_p75,
muB_p975 = \
np.percentile(fit['muB'],
[2.5, 25, 50, 75, 97.5])
muB_mean = fit['muB'].mean()
txtout.write('\nmuA:\n')
txtout.write(' mean =
{0:.3f}\n'.format(muA_mean))
txtout.write(' 95%CI = [{0:.3f},
{1:.3f}]\n'.format(muA_p025, muA_p975))
txtout.write(' Q1 = {0:.3f}, Median = {1:.3f}, Q3 = {2:.3f}\n'.format(
muA_p25, muA_p50, muA_p75))
txtout.write('\nmuB:\n')
txtout.write(' mean =
{0:.3f}\n'.format(muB_mean))
txtout.write(' 95%CI = [{0:.3f},
{1:.3f}]\n'.format(muB_p025, muB_p975))
txtout.write(' Q1 = {0:.3f}, Median = {1:.3f}, Q3 = {2:.3f}\n'.format(
muB_p25, muB_p50, muB_p75))
plt.title('平均値の事後分布', fontsize = 20)
sb.kdeplot(fit['muA'], label =
VarNameA)
sb.kdeplot(fit['muB'], label =
VarNameB)
plt.yticks([])
plt.legend(fontsize = 12)
plt.savefig(foutnm[:-4] + '_muAB.png')
plt.show()
mean_mu_d = fit['mu_d'].mean()
mu_d025P, mu_d25P, mu_d50P, mu_d75P,
mu_d975P = \
np.percentile(fit['mu_d'], [2.5, 25, 50, 75, 97.5])
txtout.write('\nmu_d:\n')
txtout.write(' mean =
{0:.3f}\n'.format(mean_mu_d))
txtout.write(' 95%CI = [{0:.3f},
{1:.3f}]\n'.format(mu_d025P, mu_d975P))
txtout.write(' Q1 = {0:.3f}, Median = {1:.3f}, Q3 = {2:.3f}\n'.format(
mu_d25P, mu_d50P, mu_d75P))
sb.kdeplot(fit['mu_d'])
plt.title('平均値の差の事後分布' +
'\n平均 =
{0:.2f} 95%CI = [{1:.2f},
{2:.2f}]'.format(
mean_mu_d, mu_d025P, mu_d975P), fontsize = 16)
plt.yticks([])
plt.xlabel('平均値({0}) - 平均値({1})'.format(VarNameA, VarNameB), fontsize = 14)
plt.savefig(foutnm[:-4] + '_mu_d.png')
plt.show()
mean_sgm_r = fit['sgm_r'].mean()
sgm_r025P, sgm_r25P, sgm_r50P,
sgm_r75P, sgm_r975P = \
np.percentile(fit['sgm_r'], [2.5, 25, 50, 75, 97.5])
txtout.write('\nsgm_r:\n')
txtout.write(' mean =
{0:.3f}\n'.format(mean_sgm_r))
txtout.write(' 95%CI = [{0:.3f},
{1:.3f}]\n'.format(sgm_r025P, sgm_r975P))
txtout.write(' Q1 = {0:.3f}, Median = {1:.3f}, Q3 = {2:.3f}\n'.format(
sgm_r25P, sgm_r50P, sgm_r75P))
plt.title(('標準偏差の比SD({0})/SD({1})の事後分布'.format(VarNameA, VarNameB)) +
'\n平均 =
{0:.3f}, 95%CI = [{1:.3f},
{2:.3f}]'.format(
mean_sgm_r, sgm_r025P, sgm_r975P), fontsize = 18)
sb.kdeplot(fit['sgm_r'])
plt.xlabel('SD({0})/SD({1})'.format(VarNameA,
VarNameB), fontsize = 14)
plt.tight_layout()
plt.savefig(foutnm[:-4] +
'_sgm_r.png')
plt.show()
sgmA_p025, sgmA_p25, sgmA_p50,
sgmA_p75, sgmA_p975 = \
np.percentile(fit['sgmA'],
[2.5, 25, 50, 75, 97.5])
sgmA_mean = fit['sgmA'].mean()
sgmB_p025, sgmB_p25, sgmB_p50,
sgmB_p75, sgmB_p975 = \
np.percentile(fit['sgmB'],
[2.5, 25, 50, 75, 97.5])
sgmB_mean = fit['sgmB'].mean()
txtout.write('\nsgmA:\n')
txtout.write(' mean =
{0:.3f}\n'.format(sgmA_mean))
txtout.write(' 95%CI = [{0:.3f},
{1:.3f}]\n'.format(sgmA_p025, sgmA_p975))
txtout.write(' Q1 = {0:.3f}, Median = {1:.3f}, Q3 = {2:.3f}\n'.format(
sgmA_p25,
sgmA_p50, sgmA_p75))
txtout.write('\nsgmB:\n')
txtout.write(' mean =
{0:.3f}\n'.format(sgmB_mean))
txtout.write(' 95%CI = [{0:.3f},
{1:.3f}]\n'.format(sgmB_p025, sgmB_p975))
txtout.write(' Q1 = {0:.3f}, Median = {1:.3f}, Q3 = {2:.3f}\n'.format(
sgmB_p25, sgmB_p50, sgmB_p75))
plt.title('標準偏差の事後分布', fontsize = 20)
sb.kdeplot(fit['sgmA'], label =
VarNameA)
sb.kdeplot(fit['sgmB'], label =
VarNameB)
plt.xlabel('標準偏差', fontsize = 14)
plt.savefig(foutnm[:-4] + '_sgmAB.png')
plt.show()
#
# サンプル値の差の事後予測
#
diff = []
for ma, mb, sgma, sgmb in
zip(fit['muA'], fit['muB'], fit['sgmA'], fit['sgmB']):
va = ss.norm.rvs(ma, sgma)
vb = ss.norm.rvs(mb, sgmb)
diff.append(va - vb)
diff = np.array(diff)
p_pos_diff = (diff > 0.0).sum() /
len(diff)
sb.kdeplot(diff)
plt.title(('サンプル値の差の事後予測 {0} -
{1}'.format(VarNameA, VarNameB)) +
'\n確率（差＞0）= {0:.3f}'.format(p_pos_diff), fontsize = 16)
plt.xlabel('{0}サンプル値 - {1}サンプル値'.format(VarNameA, VarNameB), fontsize = 14)
plt.yticks([])
plt.savefig(foutnm[:-4] +
'_postD.png')
plt.show()
txtout.close()
print('\n出力ファイル{}が保存されました。\n'.format(txtoutnm))