２つの独立なデータの正規分布によるベイズ分析
このウェブサイトで説明しているスクリプトを更新した（2025.03.22）。pystan3でArvizが使えるようになったので、MCMCサンプリングのRhat値の計算を行い、MCMCサンプリングのトレース図を出力するようにした。
Minicondaにおけるライブラリの更新は、使用する仮想環境をactivateした状態でコマンド
conda update --all
を実行すればよいが、pystanの更新はコマンド
pip install –U
pystan
を実行すればよい。アンインストールしてから、インストールをやり直してもよい。
正規分布は次式で与えられる（Gelman et
al., 2014）。
正規分布によるベイズ分析を行うスクリプトをリスト１とリスト２のように用意した。スクリプトは、リスト１とリスト２に分割した。これは、PyStan3はWindows上で直接実行できないので、WindowsにWSL(Windows Subsystem for Linux)をインストールしてWSL上で実行する必要があり、グラフィック（GUI）がWSLでは利用できなかった（2025.02）のでグラフ描画はWindows上で行ったからである。ファイル類は、まとめて圧縮ファイルupdatefiles.zipとした。
リスト１のスクリプトを、miniconda/Ubuntu24/WSL/Windows10のPython3.12仮想環境で実行した。Minicondaの使用は、Anacondaと同様である。このウェブサイトを参照されたい。
リスト１のスクリプトを実行すると、入力データファイル名の設定が求められる。
(py312) *****/updatefiles$
python ba_normal_2503.py
Input data file(*.csv) =
入力データファイルは、CSVファイルとしてスクリプトファイルと同じフォルダ内に用意する。
図１ CSVファイル（data.csv）
1行目に変数名を設定する。１列目は空欄でもよい。データは２行目から設定する。１列目はデータの識別用であり、分析では用いられない。２列目と３列目にデータ値を設定する。データ数が２つのグループで異なるときは、少ない方のデータは空欄にしておけばよい。空欄はスクリプトの実行時に無視される。
作成したデータファイルは、CSVファイルとして保存する。データファイルとリスト１およびリスト２のファイルは、同じフォルダ内に置いておく。
次の例では、図１のファイル名data.csvが設定されている。
(py312) *****/updatefiles$
python ba_normal_2503.py
Input data file(*.csv) = data.csv
ファイル名を設定すると、データが読み込まれて計算が始まる。
MCMCサンプリングが終了すると、MCMCサンプルの統計量が出力されて、スクリプトの実行終了となる。
fit.pkl was saved.
mean sd hdi_3% hdi_97% mcse_mean mcse_sd ess_bulk ess_tail r_hat
vmu1 0.404 0.076 0.268 0.557 0.001 0.002 3370.0 1983.0 1.0
vsgm1 0.613 0.111 0.422 0.825 0.002 0.002 2841.0 1858.0 1.0
vmu2 0.365 0.061 0.257 0.488 0.001 0.001 3486.0 1946.0 1.0
vsgm2 0.544 0.094 0.377 0.711 0.002 0.002 2625.0 1890.0 1.0
mu1 60.727 2.211 56.764 65.155 0.039 0.047 3370.0 1983.0 1.0
sgm1 8.896 1.602 6.123 11.962 0.033 0.029 2841.0 1858.0 1.0
mu2 56.168 2.866 51.083 61.916 0.049 0.059 3486.0 1946.0 1.0
sgm2 12.774 2.209 8.863 16.698 0.047 0.047 2625.0 1890.0 1.0
(py312) *****/updatefiles$
MCMCサンプルの分析
リスト１の実行によるMCMCサンプルのファイルを読み込んで分析を行うために、リスト２のスクリプトを実行する。WSLではグラフ描画ができなかったので、リスト２のスクリプトは直接Windows上でPython3.13.2により実行した。
リスト１の実行後、リスト２のスクリプトを実行すると、リスト１の実行により生成されたMCMCサンプルのファイルが読み込まれて、まずサンプルの分布グラフとトレース図が表示される（図２）。
図２
図２のWindowsを閉じると、次にmuの事後分布のグラフが表示される（図３）。
図３
パラメータ
と
の事後分布の中央値は60.72と56.19である。事後分布の代表値として中央値をとっている（Gelman et al., 2021）。
図３のWindowを閉じると、図４のグラフが表示される。
図４
平均値パラメータの差mu1-mu2の事後分布である。事後確率
が約90％であり、の方がmu2より大きいことがかなりの確率で言えると考えられる。
の中央値は4.49である。
図４のWindowを閉じると、図５のグラフが表示される。
図５
パラメータsgmの事後分布である。
図５を閉じると、図６のグラフが表示される。
図６
データのヒストグラムと、パラメータの点推定値（事後分布の中央値）による正規分布のグラフである。
リスト１ MCMCサンプリングのスクリプト例（ba_normal_2503.py）
import numpy
as np
import pandas as pd
import stan
import pickle
import arviz
as az
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) # DataFrame
of MCMC samples is saved
print('fit.pkl
was saved.')
i_data = az.from_pystan(posterior=fit,
posterior_model=sm)
with open('i_data.pkl',
'wb') as f:
pickle.dump(i_data, f)
# InferenceData
of MCMC samples is saved
print(az.summary(i_data))
# Ptint
summary of MCMC sample
リスト２ MCMCサンプルの描画（draw_normal_2503.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
import arviz
as az
"""
保存されたデータの読み込み
"""
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('i_data.pkl',
'rb') as f:
i_data
= pickle.load(f)
with open('fit.pkl',
'rb') as f:
pd_fit_frame
= pickle.load(f)
"""
MCMCサンプリングのトレース図
"""
az.plot_trace(i_data, var_names=['mu1', 'mu2', 'sgm1', 'sgm2'], figsize=(10,7))
plt.tight_layout()
plt.show()
"""
事後分布の中央値を点推定値とする
"""
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でデータ分析」丸善出版