２つの独立なデータの正規分布によるベイズ分析
正規分布は次式で与えられる（Gelman et al., 2014）。
正規分布によるベイズ分析を行うスクリプトをリスト１とリスト２のように用意した。スクリプトは、リスト１とリスト２に分割した。これは、PyStan3はWindows上で直接実行できないので、WindowsにWSL(Windows Subsystem for Linux)をインストールしてWSL上で実行する必要があり、グラフィック（GUI）がWSLでは利用できなかった（2025.02）のでグラフ描画はWindows上で行ったからである。ファイル類は、まとめて圧縮ファイルba_normal_files.zipとした。
なお、更新版をこのウェブサイトに用意した（2025.03.22）。更新版では、MCMCサンプリングのRhat値とサンプリングのトレース図が出力される。
リスト１のスクリプトを、miniconda/Ubuntu24/WSL/Windows10のPython3.12仮想環境で実行した。Minicondaの使用は、Anacondaと同様である。このウェブサイトを参照されたい。
リスト１のスクリプトを実行すると、入力データファイル名の設定が求められる。
(py312) *****/ba_normal_files$ python
ba_normal.py
Input data file(*.csv) =
入力データファイルは、CSVファイルとしてスクリプトファイルと同じフォルダ内に用意する。
図１ CSVファイル（data.csv）
1行目に変数名を設定する。１列目は空欄でもよい。データは２行目から設定する。１列目はデータの識別用であり、分析では用いられない。２列目と３列目にデータ値を設定する。データ数が２つのグループで異なるときは、少ない方のデータは空欄にしておけばよい。空欄はスクリプトの実行時に無視される。
作成したデータファイルは、CSVファイルとして保存する。データファイルとリスト１およびリスト２のファイルは、同じフォルダ内に置いておく。
次の例では、図１のファイル名data.csvが設定されている。
(py312) *****/ba_normal_files$ python
ba_normal.py
Input data file(*.csv) = data.csv
ファイル名を設定すると、データが読み込まれて計算が始まる。
MCMCサンプリングが終了すると、MCMCサンプルがファイルfit.pklに保存されたことが表示されて、スクリプトの実行終了となる。
fit.pkl was saved.
(py312) *****/ba_normal_files$
MCMCサンプルの分析
リスト１の実行によるMCMCサンプルのファイルを読み込んで分析を行うために、リスト２のスクリプトを実行する。WSLではグラフ描画ができなかったので、リスト２のスクリプトは直接Windows上でPython3.13.2により実行した。
リスト１の実行後、リスト２のスクリプトを実行すると、リスト１の実行により生成されたMCMCサンプルのファイルが読み込まれて、パラメータmuの事後分布のグラフが表示される（図２）。
図２
パラメータ
と
の事後分布の中央値は60.74と56.17である。事後分布の代表値として中央値をとっている（Gelman et al., 2021）。
図２のWindowを閉じると、図３のグラフが表示される。
図３
平均値パラメータの差mu1-mu2の事後分布である。事後確率
が約90％であり、の方がmu2より大きいことがかなりの確率で言えると考えられる。
の中央値は4.56である。
図３のWindowを閉じると、図４のグラフが表示される。
図４
パラメータsgmの事後分布である。
図４を閉じると、図５のグラフが表示される。
図５
データのヒストグラムと、パラメータの点推定値（事後分布の中央値）による正規分布のグラフである。
リスト１ MCMCサンプリングのスクリプト例（ba_normal.py）
import numpy as np
import pandas as pd
import stan
import pickle
inf_nm = input('Input data file(*.csv)
= ')
dframe = pd.read_csv(inf_nm)
# データの読み込み
with open('in_dframe.pkl', 'wb') as f:
pickle.dump(dframe, f)
# 読み込みデータの保存
df_keys = dframe.keys()
x_name = df_keys[1]
y_name = df_keys[2]
print(x_name, y_name)
x = dframe[x_name].values
x = x[np.isnan(x)==False] # 空セルデータの削除
print('X =\n', x)
y = dframe[y_name].values
y = y[np.isnan(y)==False] # 空セルデータの削除
print('Y =\n', y)
mean1 = np.mean(x)
sd1 = np.std(x)
print('mean1 =', mean1, ' sd1 =', sd1)
mean2 = np.mean(y)
sd2 = np.std(y)
print('mean2 =', mean2, ' sd2 =', sd2)
min1 = np.min(x)
max1 = np.max(x)
min2 = np.min(y)
max2 = np.max(y)
print(min1,' ', max1)
print(min2,' ', max2)
stan_code = """
data {
int
n1;
array[n1] real x;
int
n2;
array[n2] real y;
real
min1;
real
max1;
real
min2;
real
max2;
}
parameters {
real<lower=0, upper=1> vmu1;
real<lower=0, upper=1> vsgm1;
real<lower=0, upper=1> vmu2;
real<lower=0, upper=1> vsgm2;
}
transformed parameters {
real
mu1;
real<lower=0> sgm1;
real
mu2;
real<lower=0> sgm2;
/*
Parameters rescaled
*/
mu1
= min1 + (max1 - min1) * vmu1;
sgm1
= 0.0001 + (max1 - min1) * vsgm1 / 2.0;
mu2
= min2 + (max2 - min2) * vmu2;
sgm2
= 0.0001 + (max2 - min2) * vsgm2 / 2.0;
}
model {
vmu1
~ beta(1, 1);
vsgm1 ~ beta(1, 1);
vmu2
~ beta(1, 1);
vsgm2 ~ beta(1, 1);
x ~
normal(mu1, sgm1);
y ~
normal(mu2, sgm2);
}
"""
Data = {'n1':len(x), 'x':x,
'min1':min1, 'max1':max1,
'n2':len(y),
'y':y, 'min2':min2, 'max2':max2}
sm = stan.build(stan_code, data=Data)
fit = sm.sample()
pd_fit_frame = fit.to_frame()
with open('fit.pkl', 'wb') as f:
pickle.dump(pd_fit_frame,f) # MCMC samples saved
print('fit.pkl was saved.')
リスト２ MCMCサンプルの描画（draw_normal.py）
import numpy as np
import pandas as pd
import pickle
import matplotlib.pyplot as plt
import seaborn as sb
import scipy.stats as ss
"""
保存されたデータの読み込み
"""
with open('in_dframe.pkl', 'rb') as f:
dframe = pickle.load(f)
df_keys = dframe.keys()
x_name = df_keys[1]
y_name = df_keys[2]
print(x_name, y_name)
x = dframe[x_name].values
x = x[np.isnan(x)==False] # 空セルデータの削除
print('X =\n', x)
y = dframe[y_name].values
y = y[np.isnan(y)==False] # 空セルデータの削除
print('Y =\n', y)
mean1 = np.mean(x)
sd1 = np.std(x)
print('mean1 =', mean1, ' sd1 =', sd1)
mean2 = np.mean(y)
sd2 = np.std(y)
print('mean2 =', mean2, ' sd2 =', sd2)
min1 = np.min(x)
max1 = np.max(x)
min2 = np.min(y)
max2 = np.max(y)
print(min1,' ', max1)
print(min2,' ', max2)
"""
MCMCサンプルの読み込み
"""
with open('fit.pkl', 'rb') as f:
pd_fit_frame =
pickle.load(f)
"""
事後分布の中央値を点推定値とする
"""
mu1 = np.median(pd_fit_frame['mu1'])
sgm1 = np.median(pd_fit_frame['sgm1'])
mu2 = np.median(pd_fit_frame['mu2'])
sgm2 = np.median(pd_fit_frame['sgm2'])
"""
パラメータmu1とmu2の事後分布
"""
plt.title('Posterior distributions'
+f'\nmu1={mu1:.2f}, mu2={mu2:.2f}')
sb.kdeplot(pd_fit_frame['mu1'],
label=x_name)
sb.kdeplot(pd_fit_frame['mu2'],
label=y_name)
plt.xlabel('mu', fontsize=14)
plt.legend()
plt.show()
"""
パラメータmu1とmu2の差の事後分布と事後確率
"""
diff_mu1_mu2 = pd_fit_frame['mu1'] -
pd_fit_frame['mu2']
diff_mu_med = np.median(diff_mu1_mu2)
pos_idx = diff_mu1_mu2 > 0
neg_idx = diff_mu1_mu2 < 0
prop_pos = np.mean(pos_idx) # P(mu1-mu2>0)
prop_neg = np.mean(neg_idx) # P(mu1-mu2<0)
sb.kdeplot(diff_mu1_mu2)
plt.xlabel('mu1-mu2', fontsize=14)
plt.title(f'P(mu1>mu2)={prop_pos:.5f},
P(mu1<mu2)={prop_neg:.5f}' +
f' med.={diff_mu_med:.2f}')
plt.tight_layout()
plt.show()
"""
パラメータsgm1とsgm2の事後分布
"""
plt.title('Posterior distributions'
f'\nsgm1={sgm1:.1f},
sgm2={sgm2:.2f}')
sb.kdeplot(pd_fit_frame['sgm1'], label=x_name)
sb.kdeplot(pd_fit_frame['sgm2'],
label=y_name)
plt.xlabel('sgm', fontsize=14)
plt.legend()
plt.show()
"""
ヒストグラムと密度関数
"""
plt.hist(x, density=True, color='b',
alpha=0.3, label=x_name)
plt.hist(y, density=True, color='r',
alpha=0.3, label=y_name)
minxy = min1 if min1 < min2 else
min2
maxxy = max1 if max1 > max2 else
max2
x_coord = np.linspace(minxy, maxxy,
1000)
y_coord1 = ss.norm.pdf(x_coord,
loc=mu1, scale=sgm1)
plt.plot(x_coord, y_coord1, color='b',
label=x_name)
y_coord2 = ss.norm.pdf(x_coord,
loc=mu2, scale=sgm2)
plt.plot(x_coord, y_coord2, color='r',
label=y_name)
plt.title(f'{x_name}: mu1={mu1:.1f}, sgm1={sgm1:.1f}' +
f'\n{y_name}:
mu2={mu2:.1f},
sgm2={sgm2:.1f}')
plt.legend()
plt.show()
ベイズ分析の基本的考え方の説明を岡本（2019）、pp.169-177、岡本（2014）、pp.227-232で行っている。
Gelman, A., Carlin, J. B., Stern, H. S., Dunson, D. B., Vehtari, A., & Rubin, D. B. (2014). Bayesian data analysis, 3rd. Ed. CRC Press.
Gelman, A., Hill, J., & Vehtari, A. (2021). Regression and other stories. Cambridge University
岡本安晴（2019）「いまさら聞けないPythonでデータ分析」丸善出版