独立なサンプルの平均値の差のベイズ分析
2025.11：CmdStanPy対応に改訂
２つの独立なサンプルの分布の差のベイズ分析を行う。
ベイズ分析を行うソフトとしてStanがある。CmdStanPyはStanのPythonインターフェースであるが、CmdStanPyの簡単な説明をこのウェブサイトで行った。
データは、正規分布に従うという仮定を置いているが、分散は等しくなくてもよい。正規分布以外の分布として対数正規分布がよく知られているが、対数正規分布の場合のベイズ分析の例として、２パラメータモデルの場合、３パラメータモデルの場合を用意した。
分析例のデータとして、A市から１５人、B市から２０人をランダムに選び血圧を測定した結果（仮想データ）である図１のものを用いる。
図１ data.xlsx
データの形式は、１行目に変数名を書き、２行目以降にデータ値が書かれている。第１列は各データの識別用で数値でも文字列でもよい。２列目に第１グループ（図１の場合はA市）のデータを書き、３列目に第２グループ（図１の場合はB市）のデータを書く。人数が揃わないときは、少ない方のデータは欠損値として空欄にしておく。
図１の形式のデータに対して以下の確率モデルを設定して分布の差のベイズ分析を行う。
上のモデルに対するStanスクリプト例をリスト１としてこのウェブサイトの後半に示す。なお、サンプルファイルは圧縮ファイルfilesBysDifMuInd.zipとしてまとめた。ダウンロード解凍すれば自由に利用できる。
リスト１のStanスクリプトを利用するPythonスクリプトをリスト２としてこのウェブサイトの後半に示す。
リスト１のStanスクリプトファイル（IndDiffMu2Sgm.stan）、リスト２のPythonスクリプトファイル（DiffMeansInd.py）、および入力データファイル」（data.xlsx）を同じフォルダに置く。そのフォルダにカレントディレクトリ（フォルダ）を移し、CmdStanPyのインストールされた環境で次のコマンドを実行する。CmdStanPyの簡単な説明を、このウェブサイトで行った。
(stan) *****/filesBysDifMuInd$ python DiffMeansInd.py
上のコマンドを実行すると、入力データファイル名の設定を求めるダイアログウィンドウが表示される（図２）。
図２
入力データファイルを設定する。図２では、図１のファイル名data.xlsxが設定されている。データファイル名を設定して「Open」ボタンをクリックすると、計算結果などが書き出される出力ファイル名の設定を求めるダイアログウィンドウが表示される（図３）。
図３
適当なテキストファイル名を設定する。ファイル名を設定して「Save」ボタンをクリックすると、次に、グラフ画像を保存するファイルの名前の設定を求めるダイアログウィンドウが表示される（図４）。
図４
画像は複数枚保存されるので、図４のダイアログウィンドウで設定した名前を基に複数の名前が作成されて保存される。設定した名前が図４のように「MyFig.png」であれば、「MyFig_XXX.png」という形式の名前（XXXの部分は図によって変わる）で画像が保存される。保存ファイル名を設定して「Save」ボタンをクリックすると、Stanスクリプトのコンパイルが始まる。コンパイルは多少時間が掛かる。コンパイル後、MCMCサンプリングが始まる。サンプリングが終了すると、トレース図が表示される（図５）。
図５
左側に事後分布のグラフ、右側にトレース図が表示されている。
図５のWindowを閉じると、図６のグラフが表示される。
図６
平均値パメータと標準偏差パラメータの事後分布である。
左上のグラフは、2つのデータそれぞれの平均値パラメータの事後分布を表し、A市の分布の方がはっきりと左側にある。
右上のグラフは、「A市の平均―B市の平均」の事後分布のグラフであり、95％CIが負の領域にある。
これらの結果は、A市の平均の方がB市の平均より低いことを示すものである。
左下のグラフは、A市およびB氏の標準偏差の事後分布のグラフである。ほぼ重なっている。
右下のグラフは、比「A市の標準偏差／B市の標準偏差」の事後分布のグラフである。95％CIの中ほどに１があり、比が１であることが強く支持されている。すなわち、A市の標準偏差とB市の標準偏差は同じである（等分散）と考えてよさそうである。
図６のWindowを閉じると、図７の事後予測分布のグラフが表示される。
図７
２つの群（A市とB市）からそれぞれランダムにサンプルを採り、その差を求めたものである。２つの群は、パラメータ値が事後分布から取られた値である正規分布に従うものとしてシミュレーションによって求めている。スクリプトは次のようである。
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)
A市とB市からランダムに１人ずつ選んで比べたとき、A市の人の方が値が高いのは100回中数回のレベルであるという事後予測である。
図７のWindowを閉じると、プログラムの実行終了である。
実行終了後、図３で指定した出力ファイルを開くと以下のようになっている。PythonのIDLEでは開けないことがあるが、そのときはメモ帳やワードパッドなど他のエディタで開けばよい。
入力データファイル:
/mnt/c/WorkSpace/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]
Mean
MCSE StdDev MAD ...
ESS_bulk
ESS_tail
ESS_bulk/s R_hat
lp__ -62.291300 0.038533 1.510790 1.263430 ... 1695.80 2390.06 32611.5 1.00163
muA 125.042000 0.019111 1.060910 0.962571 ... 3226.49 2353.34 62048.0 1.00255
muB 133.944000 0.014595 0.911423 0.859656 ... 3955.50 2656.55 76067.3 1.00075
sgmA 3.938680 0.015041 0.849723 0.755009 ... 3816.47 2319.37 73393.7 1.00051
sgmB 4.176820 0.012027 0.730190 0.674013 ... 4045.82 2807.16 77804.2 1.00029
mu_d -8.901650 0.024710 1.397440 1.313070 ... 3278.92 2529.74 63056.2 1.00174
sgm_r 0.971908 0.004631 0.274173 0.242892 ... 3747.24 2751.21 72062.2 1.00084
[7 rows x 11 columns]
muA:
mean = 125.042
95%CI = [122.852, 127.107]
Q1 = 124.398, Median = 125.060, Q3 = 125.696
muB:
mean = 133.944
95%CI = [132.144, 135.780]
Q1 = 133.356, Median = 133.946, Q3 = 134.518
mu_d:
mean = -8.902
95%CI = [-11.689, -6.192]
Q1 = -9.780, Median = -8.877, Q3 = -8.012
sgm_r:
mean = 0.972
95%CI = [0.560, 1.621]
Q1 = 0.782, Median = 0.927, Q3 = 1.117
sgmA:
mean = 3.939
95%CI = [2.696, 6.072]
Q1 = 3.339, Median = 3.823, Q3 = 4.369
sgmB:
mean = 4.177
95%CI = [3.000, 5.900]
Q1 = 3.659, Median = 4.083, Q3 = 4.592
リスト１ Stanスクリプト（IndDiffMu2Sgm.stan）
data {
int NA;
int NB;
array[NA] real DataA;
array[NB] real DataB;
}
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スクリプト（DiffMeansInd.py）
import tkinter
as tk
import tkinter.filedialog
as fd
import pandas as pd
import numpy
as np
from cmdstanpy
import CmdStanModel
import arviz
as az
import matplotlib.pyplot
as plt
import japanize_matplotlib
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('Input
Data File', types = [('Excel(*.xlsx)', '*.xlsx')])
txtoutnm = askoutfilename(msg
= 'Output Text File',
types = [('Text (*.txt)', '*.txt')])
txtout = open(txtoutnm,
'w')
foutnm = askoutfilename(msg
= 'Output Image File',
types = [('Image(*.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}
model = CmdStanModel(stan_file='IndDiffMu2Sgm.stan')
fit = model.sample(data=StanData)
print(fit.summary())
txtout.write(f'{fit.summary()}')
infdata = az.from_cmdstanpy(fit) # Arviz InfereceData型へ変換
az.plot_trace(infdata, var_names=['muA', 'muB', 'sgmA', 'sgmB'])
plt.tight_layout()
plt.savefig(foutnm[:-4]
+ '_trace.png')
plt.show()
fit = fit.draws_pd()
# Pandas DataFrame型への変換
print(fit.keys())
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.figure(figsize=(10,7))
plt.subplot(221)
plt.title('平均値の事後分布', fontsize = 20)
sb.kdeplot(fit['muA'],
label = VarNameA)
sb.kdeplot(fit['muB'],
label = VarNameB)
plt.xlabel('mu', fontsize=14)
plt.yticks([])
plt.legend(fontsize =
12)
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))
plt.subplot(222)
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 = 14)
plt.yticks([])
plt.xlabel('平均値({0}) - 平均値({1})'.format(VarNameA, VarNameB), fontsize = 14)
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))
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.subplot(223)
plt.title('標準偏差の事後分布', fontsize =
20)
sb.kdeplot(fit['sgmA'],
label = VarNameA)
sb.kdeplot(fit['sgmB'],
label = VarNameB)
plt.xlabel('標準偏差', fontsize = 14)
plt.legend(fontsize=14)
plt.subplot(224)
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 = 14)
sb.kdeplot(fit['sgm_r'])
plt.xlabel('SD({0})/SD({1})'.format(VarNameA, VarNameB), fontsize = 14)
plt.tight_layout()
plt.savefig(foutnm[:-4]
+ 'mu_sgm.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))