２つの独立なデータの平均値の分析
ｔ分布に基づくベイズ分析
正規分布の観点から見ると外れ値と見做される値も、ｔ分布のように裾の厚い分布においては他の値と同様に扱われる（Martin, 2018）。
ｔ分布は、次式で与えられる（Gelman et al., 2014）。
は、自由度（degrees of freedom, df）と呼ばれることもあるパラメータである。
パラメータ
の幾つかの値に対するｔ分布と正規分布のグラフを図１に示す。
図１
パラメータの値が小さいとき裾の厚い分布であるが、の値が大きくなると正規分布に近付いている。
確率分布としてｔ分布を用いるベイズ分析のスクリプトを、リスト１とリスト２に分けて用意した。スクリプトファイルなどは圧縮ファイルtest_t_distri_files.zipとしてまとめた。
リスト１のスクリプトは、Windows上では直接実行できない。Windowsで実行するときは、WindowsにWSL（Windows Subsystem for Linux）をインストールして、WSL上にインストールされたUbuntu上でリスト１のスクリプトを実行する。Ubuntu上での実行は、miniconda上に用意したPython3.12仮想環境で行った。すなわち、リスト１のスクリプトは、Python3.12/miniconda/Ubuntu/WSL/Windows 10で実行した。Minicondaの使用は、Anacondaと同様である。このウェブサイトを参照されたい。
リスト１のスクリプトでMCMCサンプリングを行った後、リスト２のスクリプトでMCMCサンプルの分析を行った。リスト２のスクリプトは直接Windows上で実行した。WSL/Windowsでは描画ソフトが使えないようである（2025.02）。
リスト１のスクリプトを実行すると、入力データファイル名の設定が求められる。
(py312) yasuharu@PSY-PM-PC/mnt/*****/test_t_distri_files$
python ba_t_distri_means.py
Input data file(*.csv) =
入力データファイルは、図２の形式のCSVファイルとして用意する。
・
・
・
図２ 入力データファイル例sampleN100.csv
１行目に変数名を書き、２行目からデータ値を書く。１列目はデータの識別用で、プログラムでは用いられない。２列目と３列目に分析データ値を書くが、データ数が異なるときは少ない方のデータのセルは空になる（図２）。
次の例では、図２のファイル名が設定されている。
(py312) yasuharu@PSY-PM-PC:/mnt/*****/test_t_distri_files$
python ba_t_distri_means.py
Input data file(*.csv) =
sampleN100.csv
ファイル名を設定すると、データが読み込まれて計算が始まる。
MCMCサンプリングが終了すると、MCMCサンプルがファイルに保存されて、スクリプトの実行終了となる。
fit.pkl was saved.
(py312) yasuharu@PSY-PM-PC:/mnt/*****/ba_t_distri/test_t_distri_files$
リスト１のスクリプトの実行後、リスト２のスクリプトを実行すると、MCMCサンプルのファイルを読み込んで分析が行われる。
まず、パラメータmu1とmu2の事後分布のグラフが表示される（図３）。
図３
分布の代表値は中央値である（Gelman et al., 2021）。
図３のWindowを閉じると、パラメータの差mu1-mu2の事後分布のグラフが表示される（図４）。
図４
タイトル部にmu1>mu2の事後確率が1.000であること、mu1<mu2の事後確率が0.000であることが表示されている。ほぼ確実にmu1の方がmu2より大きいと考えられる。
図４のWindowを閉じると、パラメータsgmの事後分布が表示される。
図５
図５のWindowを閉じると、パラメータnuの事後分布のグラフが表示される（図６）。
図６
KDEによる推定分布では負の領域にも分布が正の値をとっている領域があるが、ヒストグラムではnuの分布は正の領域に限られていることが分かる。
図６のWindowを閉じると、図７のグラフが表示される。
図７
データのヒストグラムと、事後分布の中央値をパラメータとするｔ分布のグラフが描かれている。
リスト１ ベイズ分析のＭＣＭＣサンプリングを行うスクリプト（ba_t_distri_means.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;
real<lower=0, upper=1> vnu1;
real<lower=0, upper=1> vnu2;
}
transformed parameters {
real
mu1;
real<lower=0> sgm1;
real
mu2;
real<lower=0> sgm2;
real
nu1;
real
nu2;
/*
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;
nu1
= 0.0001 + vnu1*100;
nu2
= 0.0001 + vnu2*100;
}
model {
vmu1
~ beta(1, 1);
vsgm1 ~ beta(1, 1);
vmu2
~ beta(1, 1);
vsgm2 ~ beta(1, 1);
vnu1
~ beta(1, 1);
vnu2 ~ beta(1, 1);
x ~ student_t(nu1, mu1, sgm1);
y ~ student_t(nu2, 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_t_distri_means.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)
#max_v = np.max(x) if np.max(x) > np.max(y) else np.max(y);
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'])
nu1 = np.median(pd_fit_frame['nu1'])
mu2 = np.median(pd_fit_frame['mu2'])
sgm2 = np.median(pd_fit_frame['sgm2'])
nu2 = np.median(pd_fit_frame['nu2'])
"""
パラメータ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}')
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()
"""
パラメータnu1とnu2の事後分布
"""
plt.title('Posterior distributions' +
f'\nnu1={nu1:.2f},
nu2={nu2:.2f}')
plt.hist(pd_fit_frame['nu1'],
color='g', alpha=0.5,
density=True, label=x_name)
plt.hist(pd_fit_frame['nu2'],
color='orange', alpha=0.5,
density=True, label=y_name)
sb.kdeplot(pd_fit_frame['nu1'],
color='g', label=x_name)
sb.kdeplot(pd_fit_frame['nu2'],
color='orange', label=y_name)
plt.xlabel('nu', fontsize=14)
plt.legend()
plt.show()
"""
データのヒストグラムとｔ分布のグラフ
"""
plt.hist(x, density=True, color='b',
alpha=0.5, label=x_name)
plt.hist(y, density=True, color='orange',
alpha=0.5, label=y_name)
"""
グラフの横軸の範囲
"""
min_xy = min1 if min1 < min2 else min2
max_xy = max1 if max1 > max2 else max2
coord_xy = np.linspace(min_xy, max_xy, 1000)
density1 = ss.t.pdf(coord_xy, nu1, loc=mu1,
scale=sgm1) # t分布
density2 = ss.t.pdf(coord_xy, nu2, loc=mu2,
scale=sgm2) # t分布
plt.plot(coord_xy,
density1, c='g', label=x_name)
plt.plot(coord_xy,
density2, c='r', label=y_name)
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.
Martin, O. (2018). Bayesian analysis with Python, 2nd
Ed. Packt
岡本安晴（2019）「いまさら聞けないPythonでデータ分析」丸善出版