２つの独立なデータの正規分布によるベイズ分析
CmdStanPy対応：2026.03改訂
正規分布は次式で与えられる（Gelman et al., 2014）。
独立な２のデータ、XとY、が与えられたとき、それらが正規分布に従うとしてベイズ分析を行うStanスクリプトをリスト１のように用意した、パラメータは、データXに関するものは添字１を付け、データYに関するものは添字２を付けている。ファイルはbanormalcspfiles.zipにまとめた。
リスト１ 正規分布に従う２つの独立なデータのStanスクリプト（ba_normal.stan）
data {
int n1;
array[n1] real x;
int n2;
array[n2] real y;
}
transformed data {
real min1;
real min2;
real max1;
real max2;
min1 = min(x);
min2 = min(y);
max1 = max(x);
max2 = max(y)
;}
parameters {
real<lower=min1,
upper=max1> mu1;
real<lower=0.0001,
upper=max1-min1> sgm1;
real<lower=min2,
upper=max2> mu2;
real<lower=0.0001,
upper=max2-min2> sgm2;
}
transformed parameters {
}
model {
mu1 ~ uniform(min1, max2);
sgm1 ~ uniform(0.0001,
max1-min1);
mu2 ~ uniform(min2, max2);
sgm2 ~ uniform(0.0001,
max2-min2);
x ~ normal(mu1, sgm1);
y ~ normal(mu2, sgm2);
}
事前分布は、十分に広い範囲で一様分布としている。すなわち、
mu1 ~ uniform(min1, max2);
sgm1 ~ uniform(0.0001,
max1-min1);
mu2 ~ uniform(min2, max2);
sgm2 ~ uniform(0.0001,
max2-min2);
である。Stanは随分使い易くなったと思う（cmdstan 2.38.0）。
リスト１のStanスクリプトをCmdStanPyを利用して用い、ベイズ分析を行うPythonスクリプトをリスト２のように用意した。CmdStanPyの簡単な準備方法の説明を、このウェブサイトで行った。
リスト２ リスト１のStanスクリプトをCmdStanPyを利用して用いるPythonスクリプト（ba_normal.py）
from cmdstanpy
import CmdStanModel
import numpy
as np
import scipy.stats
as ss
import matplotlib.pyplot
as plt
import arviz
as az
import pandas as pd
import seaborn as sb
inf_nm = input('Input data file(*.xlsx) = ')
dframe = pd.read_excel(inf_nm)
# データの読み込み
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)
Data = {'n1':len(x), 'x':x, # 'min1':min1, 'max1':max1,
'n2':len(y), 'y':y} #, 'min2':min2, 'max2':max2}
model = CmdStanModel(stan_file="ba_normal.stan") # コンパイル
fit = model.sample(data=Data)
# サンプリング
print(fit.summary())
inf_data = az.from_cmdstanpy(fit)
# Arviz(az)用のデータに変換
# トレースプロット
az.plot_trace(inf_data, var_names=['mu1','sgm1','mu2','sgm2'])
plt.tight_layout()
plt.show()
pdf_fit = fit.draws_pd()
# Pandas DataFrameに変換
"""
事後分布の中央値を点推定値とする
"""
mu1 = np.median(pdf_fit['mu1'])
sgm1 = np.median(pdf_fit['sgm1'])
mu2 = np.median(pdf_fit['mu2'])
sgm2 = np.median(pdf_fit['sgm2'])
"""
パラメータmu1とmu2の事後分布
"""
plt.title('Posterior distributions'
+f'\nmu1={mu1:.2f}, mu2={mu2:.2f}')
sb.kdeplot(pdf_fit['mu1'],
label=x_name)
sb.kdeplot(pdf_fit['mu2'],
label=y_name)
plt.xlabel('mu', fontsize=14)
plt.legend()
plt.show()
"""
パラメータmu1とmu2の差の事後分布と事後確率
"""
diff_mu1_mu2 = pdf_fit['mu1']
- pdf_fit['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(pdf_fit['sgm1'],
label=x_name)
sb.kdeplot(pdf_fit['sgm2'],
label=y_name)
plt.xlabel('sgm', fontsize=14)
plt.legend()
plt.show()
"""
パラメータsgm1とsgm2の差の事後分布と事後確率
"""
diff_sgm1_sgm2 = pdf_fit['sgm1']
- pdf_fit['sgm2']
diff_sgm_med = np.median(diff_sgm1_sgm2)
prop_pos = np.mean(diff_sgm1_sgm2
> 0) # P(mu1-mu2>0)
prop_neg = np.mean(diff_sgm1_sgm2
< 0) # P(mu1-mu2<0)
sb.kdeplot(diff_sgm1_sgm2)
plt.xlabel('sgm1-sgm2', fontsize=14)
plt.title(f'P(sgm1>sgm2)={prop_pos:.5f},
P(sgm1<sgm2)={prop_neg:.5f}' +
f' med.={diff_sgm_med:.2f}')
plt.tight_layout()
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()
リスト１とリスト２のスクリプトファイル、および入力データファイルを同じフォルダに置き、CmdStanPyのインストールされた環境においてリスト２のスクリプトを実行すると、入力データファイル名の設定が求められる。
(stan) ****/banormalcspfiles$
python ba_normal.py
Input data file(*.xlsx) = data.xlsx
上の例では、data.xlsxが設定されている。入力データファイルは、図１a、図１ｂの形式で、Excelファイルとして用意する。
図１a
図１ｂ
ファイルの形式は、1行目に変数名を設定する。１列目は空欄でもよい。データは２行目から設定する。１列目はデータの識別用であり、分析では用いられない。２列目と３列目にデータ値を設定する。データ数が２つのグループで異なるときは、少ない方のデータは空欄にしておけばよい。空欄は、スクリプトの実行時に次のコードにより削除される。
x = x[np.isnan(x)==False] # 空セルデータの削除
y = y[np.isnan(y)==False] # 空セルデータの削除
入力データファイル名を設定してEnterキーを押すと、データが読み込まれ、Stanスクリプトのコンパイルが始まる。コンパイルには、時間が掛かるが、コンパイル後、ＭＣＭＣサンプリングが実行される。
MCMCサンプリングが終了すると、トレース図が表示される（図２）。
図２
図２のWindowの右上角のX印アイコンをクリックして閉じると、次に、パラメータmu1とmu2の事後分布が表示される（図３）。
図３
図３のWindowを閉じると、図４のグラフが表示される（図４）。
図４
パメータの差mu1-mu2の事後分布のグラフである。mu1-mu2の正である確率が89％であるので、かなりの確信をもってmu1-mu2が正であると言える。
図４のWindowを閉じると、図５のグラフが表示される。
図５
sgm1とsgm2の事後分布のグラフである。
図５のWindowを閉じると、図６のグラフが表示される。
図６
パラメータの差sgm1-sgm2の事後分布と、差sgm1-sgm2の事後分布が正である確率および負である確率が表示されている。正である確率が7.55％であるので、sgm1の方が小さい傾向にあると言える。
図６のWindowを閉じると、図７のグラフが表示される。
図７
データのヒストグラムと、パラメータの点推定値によるモデルのグラフである。点推定値は、中央値を採用している。Gelman,ら(2021)は、平均値より中央値の方が安定しているとして、中央値の採用を勧めている。
図７のWindowを閉じると、スクリプトの実行終了である。端末には、以下のような出力が表示されている。
14:59:17 - cmdstanpy
- INFO - CmdStan done processing.
Mean
MCSE StdDev MAD
5%
50%
95% ESS_bulk ESS_tail ESS_bulk/s R_hat
lp__ -96.99250 0.042301 1.59494 1.38308 -100.14300 -96.63370
-95.1223 1603.31 2098.07 36438.8 1.00024
mu1 60.69090 0.036420 2.17851 2.02340 57.12600 60.67320 64.4135 3689.83 2105.53 83859.8 1.00029
sgm1 8.98809 0.037310 1.82228 1.66992 6.55722 8.70664 12.3126 3350.80 1785.96 76154.6 1.00006
mu2 56.11620 0.053874 3.01262 2.91021 51.19650 56.13930 60.9896 3253.98 2081.17 73954.0 1.00144
sgm2 12.89950 0.043120 2.27494 2.17508 9.72999 12.57620 17.0879 3209.08 2300.47 72933.7 1.00027
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