ベイズ分析における点推定値例のRMSEとbiasによる比較
――実験心理学の観点からの検討――
基礎心理学研究、第41巻、第１号、1-7
https://doi.org/10.14947/psychono.41.1
要約
ベイズ分析におけるパラメータの点推定値としての事後分布の平均値，中央値およびMAP推定値を，真値との関係を表すbiasとRMSE（２乗平均平方根誤差）の２つの指標で比較した。全体としての傾向は，MAP推定値が他の推定値よりも真値に近いという結果であった。比較は，２項分布の比率パラメータ，1変量正規分布の標準偏差パラメータ，２変量正規分布の相関係数パラメータにおいて行われた。これらの比較は単純なモデルについてであるが，MAP推定値は任意の最高密度区間に含まれていること，およびフラットな事前分布に対するMAP推定値は漸近的に真値に収束する最尤推定値に一致することを考えれば，MAP推定値は一般的に真値の推定値として適していることが示唆される。MCMCサンプリングからMAP推定値を算出する簡単なPythonスクリプトも示した。
本文中の図および表の作成に用いたプログラムを公開用に整理したものを以下に示す。これらのファイルは圧縮ファイルPointEstimates.zipにまとめたので、自由に利用されたい。
Figure 1のベータ分布の描画プログラム
リスト１のプログラムを実行すると、Figure 1のグラフを得る。
Figure 1
Figures 2, 3 and 4 のグラフ描画プログラム
リスト２のプログラムを実行すると、Figures 2, 3, 4を得る。
Figure 2
Figure 3
Figure 4
Table 1のデータの出力プログラム
リスト３のプログラムを実行すると以下の出力ファイルResultsCheck.txtを得る
*/n_samples
20
*/Biases_mu_mean
-0.007041961545885269
*/Biases_mu_median
-0.006976374519543286
*/Biases_mu_map
-0.007209074323586536
*/SEBiases_mu_mean
0.007150514505840424
*/SEBiases_mu_median
0.0071495556400130765
*/SEBiases_mu_map
0.0071892461915072006
*/Biases_sgm_mean
0.05242901233745667
*/Biases_sgm_median
0.027014065601938864
*/Biases_sgm_map
-0.014603021898702995
*/SEBiases_sgm_mean
0.0053799409384119945
*/SEBiases_sgm_median
0.005250154481094522
*/SEBiases_sgm_map
0.005062170454121036
*/rms_mu_mean
0.22611571609012457
*/rms_mu_median
0.22608339060008384
*/rms_mu_map
0.22734455413293206
*/se_rms_mu_mean
0.04788157166709367
*/se_rms_mu_median
0.04786410713297506
*/se_rms_mu_map
0.04789616874476579
*/rms_sgm_mean
1.066077692640595
*/rms_sgm_median
1.0403338160886175
*/rms_sgm_map
0.9983021329177043
(*/se_rms_sgm_mean
0.1074652838233129
*/se_rms_sgm_median
0.10488630620335115
*/se_rms_sgm_map
0.10089430461605017
これらの値は、乱数を用いたシミュレーションによるものなので、プログラムの実行ごとに異なる値が得られるが、SEの推定値からランダムな変動の大きさが推測できる。
Tables 2 and 3のデータの出力プログラム
リスト５のプログラムを実行すると出力ファイルResultsN50.txtを得る。
*/bias_mean
0.009920408994082704 0.03211452354398249
0.028884218740121793 -0.0026328761413996504 -0.026355866314786165
-0.012562849096223934 -0.011020725687830983
*/SE_bias_mean
0.002166293613514624 0.006352606295812506
0.010099708755636351 0.009203343499895711 0.009020196954548528
0.006314623678249787 0.0022096789585528158
*/bias_median
0.006090627925731067 0.024341765670074733
0.024130178778043306 -0.0028000087042905856 -0.021459953799896535
-0.004869294532790663 -0.007182361311745917
*/SE_bias_median
0.0021051499843991453 0.00635456231224322
0.010251353698944778 0.00938546793307021 0.00915324564903438
0.006312780165466433 0.0021495696801992837
*/bias_MAP
-0.0009005450272513482
0.008310915545777241 0.014341217060853001 -0.0029016450822541505
-0.011427571378569002 0.010712035601780039 -0.0003035151757588828
*/SE_bias_MAP
0.0020021970253489233
0.0063518630703183025 0.01063789456229552 0.009772134749041317
0.009408831241654925 0.0062622276474516075 0.002017314661103084
*/RMSE_mean
0.032129227979843716 0.09519510683657403
0.14537234044003708 0.12985583090725422 0.1299463672284583 0.0899602420227911
0.03306221307717857
*/SE_RMSE_mean
0.011974161331843314 0.030981311084084282
0.044360908885495584 0.04001050562758937 0.04044402305213406
0.028420072573137268 0.012871139613899383
*/RMSE_median
0.03031493663915568 0.09288828521289981
0.14661250222055833 0.1324279226785622 0.13089357850366642 0.08918574728806637
0.031162407259530567
*/SE_RMSE_median
0.011252491414637482 0.030074829818920106
0.04466027762934707 0.04078827914170566 0.040889195247272775
0.028106586876200927 0.012099536376119631
*/RMSE_MAP
0.028258817634165195 0.08998865509362647
0.1507496794350532 0.1378834596219362 0.1332189338985517 0.08898669132699984 0.02845934383291899
*/SE_RMSE_MAP
0.010333327247190044 0.028340088036870655
0.04603693301893274 0.04251454110788044 0.04156855256467426
0.027453424591980552 0.01051948475436297
各行は、相関係数パラメータの値
が−0.9から＋0.9の範囲に対する結果が並んで出力されている。乱数を用いたシミュレーション結果であるので、プログラムの実行ごとに結果はランダムに変動するが、SEの推定値から変動の大きさが推測できる。
import numpy as np
import scipy.stats as ss
import matplotlib.pyplot as plt
import scipy.optimize as so
class Func_med:
def __init__(self, f):
self.f = f
def v(self, t):
#return self.f(t[0]) - 0.5
return self.f(t) - 0.5
X = 1
N = 100
d_func = ss.beta(X+1, N-X+1)
t_range = np.linspace(0, 0.05)
t_v = d_func.pdf(t_range)
mean_v = (X+1) / (N+2)
mode_v = X / N
print('mean = ', mean_v)
print('mode = ', mode_v)
d_func = ss.beta(X+1, N-X+1)
func_med = Func_med(d_func.cdf)
f_med = func_med.v
median_v = so.bisect(f_med, 0.0, 1.0)
plt.title('Beta(θ|2,100)', fontsize = 16)
plt.plot([0, 0.05], [0.0, 0.0], color
= 'k', ls = '-', lw = 1)
plt.plot(t_range, t_v, color = 'k', lw
= 2)
plt.plot([mean_v, mean_v], [0,
d_func.pdf(mean_v)], c = 'k', lw = 1.5,
ls = '-.', label = f'mean:{mean_v:.4f}')
plt.plot([median_v, median_v], [0,
d_func.pdf(median_v)], c = 'k', lw = 1.5,
linestyle = (0, (5,1,2,1,2,1)), label = f'median:{median_v:.4f}')
plt.plot([mode_v, mode_v], [0,
d_func.pdf(mode_v)], c = 'k',
ls = '--', lw = 1.5, label = f'mode:{mode_v:.4f}')
plt.xlabel('θ', fontsize = 16)
plt.xticks(fontsize = 14)
plt.yticks([])
plt.legend(fontsize = 14)
plt.tight_layout()
plt.savefig('Figure1.png')
plt.show()
import numpy as np
import scipy.stats as ss
import matplotlib.pyplot as plt
import scipy.optimize as so
X = 1
N = 100
d_func = ss.beta(X+1, N-X+1)
t_range = np.linspace(0, 0.05)
t_v = d_func.pdf(t_range)
mean_v = (X+1) / (N+2)
mode_v = X / N
print('mean = ', mean_v)
print('mode = ', mode_v)
class Func_med:
def __init__(self, f):
self.f = f
def v(self, t):
return self.f(t) - 0.5
N = 100
p_true_range = [v/100 for v in
range(1, 100)]
print(p_true_range)
bias_mean = []
bias_median = []
bias_mode = []
RMSE_mean = []
RMSE_median = []
RMSE_mode = []
RMSER_mean = []
RMSER_median = []
RMSER_mode = []
for p_true in p_true_range:
print(p_true, end = '\r')
v_sum_mean = 0.0
v_sum_mode = 0.0
v_sum_med = 0.0
ssum_mean = 0.0
ssum_mode = 0.0
ssum_med = 0.0
for X in range(0, N + 1):
p =
ss.binom.pmf(X, N, p_true)
#print('p =', p)
v_sum_mean += p * (X + 1) / (N + 2)
v_sum_mode += p * X / N
ssum_mean += p * ((X +1) / (N+2) - p_true) ** 2
ssum_mode += p * ((X / N) - p_true) ** 2
d_func = ss.beta(X+1, N-X+1)
func_med
= Func_med(d_func.cdf)
f_med = func_med.v
v_med = so.bisect(f_med, 0.0, 1.0)
v_sum_med += p * v_med
ssum_med += p * (v_med - p_true) ** 2
bias_mean.append(v_sum_mean
- p_true)
bias_median.append(v_sum_med
- p_true)
bias_mode.append(v_sum_mode
- p_true)
RMSE_mean.append(ssum_mean
** 0.5)
RMSE_median.append(ssum_med
** 0.5)
RMSE_mode.append(ssum_mode
** 0.5)
RMSER_mean.append((ssum_mean
** 0.5) / p_true)
RMSER_median.append((ssum_med
** 0.5) / p_true)
RMSER_mode.append((ssum_mode
** 0.5) / p_true)
N2 = 200
p_true_range = [v/100 for v in
range(1, 100)]
print(p_true_range)
bias_mean2 = []
bias_median2 = []
bias_mode2 = []
RMSE_mean2 = []
RMSE_median2 = []
RMSE_mode2 = []
RMSER_mean2 = []
RMSER_median2 = []
RMSER_mode2 = []
for p_true in p_true_range:
print(f'{p_true} for N2',
end = '\r')
v_sum_mean = 0.0
v_sum_mode = 0.0
v_sum_med = 0.0
ssum_mean = 0.0
ssum_mode = 0.0
ssum_med = 0.0
for X in range(0, N2 + 1):
p =
ss.binom.pmf(X, N2, p_true)
v_sum_mean += p * (X + 1) / (N2 + 2)
v_sum_mode += p * X / N2
ssum_mean += p * ((X +1) / (N2+2) - p_true) ** 2
ssum_mode += p * ((X / N2) - p_true) ** 2
d_func =
ss.beta(X+1, N2-X+1)
func_med = Func_med(d_func.cdf)
f_med = func_med.v
v_med = so.bisect(f_med, 0.0, 1.0)
v_sum_med += p * v_med
ssum_med += p * (v_med - p_true) ** 2
bias_mean2.append(v_sum_mean
- p_true)
bias_median2.append(v_sum_med - p_true)
bias_mode2.append(v_sum_mode
- p_true)
RMSE_mean2.append(ssum_mean
** 0.5)
RMSE_median2.append(ssum_med
** 0.5)
RMSE_mode2.append(ssum_mode
** 0.5)
RMSER_mean2.append((ssum_mean
** 0.5) / p_true)
RMSER_median2.append((ssum_med ** 0.5) / p_true)
RMSER_mode2.append((ssum_mode ** 0.5) / p_true)
plt.title('Bias', fontsize = 16)
plt.plot(p_true_range, bias_mean, c =
'k', ls = (0,(4,1)), lw = 2, label = 'mean(N=100)')
plt.plot(p_true_range, bias_median, c
= 'k', ls = (0,(4,1,2,1)), lw = 2, label = 'median(N=100)')
plt.plot(p_true_range, bias_mean2, c =
'k', ls = (0,(4,1)), lw = 1, label = 'mean(N=200)')
plt.plot(p_true_range, bias_median2, c
= 'k', ls = (0,(4,1,2,1)), lw = 1, label = 'median(N=200)')
plt.plot(p_true_range, bias_mode, c =
'k', ls = (0,(1,2)), lw = 2, label = 'MAP(N=100)')
plt.plot(p_true_range, bias_mode2, c =
'k', ls = (2,(2,2)), lw = 1, label = 'MAP(N=200)')
plt.xticks(fontsize = 14)
plt.xlabel(r'$\theta_0$', fontsize =
16)
plt.legend()
plt.tight_layout()
plt.savefig('Figure_2.png')
plt.show()
plt.title('RMS Error', fontsize = 16)
plt.plot(p_true_range, RMSE_mode, c =
'k', ls = (0,(1,2)), lw = 2, label = 'MAP(N=100)')
plt.plot(p_true_range, RMSE_median, c
= 'k', ls = (0,(4,1,2,1)), lw = 2, label = 'median(N=100)')
plt.plot(p_true_range, RMSE_mean, c =
'k', ls = (0,(4,1)), lw = 2, label = 'mean(N=100)')
plt.plot(p_true_range, RMSE_mode2, c =
'k', ls = (2,(2,2)), lw = 1, label = 'MAP(N=200)')
plt.plot(p_true_range, RMSE_median2, c
= 'k', ls = (0,(4,1,2,1)), lw = 1, label = 'median(N=200)')
plt.plot(p_true_range, RMSE_mean2, c =
'k', ls = (0,(4,1)), lw = 1, label = 'mean(N=200)')
plt.legend(fontsize = 12)
plt.xlabel(r'$θ_0$', fontsize = 13)
plt.tight_layout()
plt.savefig('Figure_3.png')
plt.show()
plt.title(r'RMS Error / $θ_0$', fontsize = 16)
plt.plot(p_true_range[:20],
RMSER_mean[:20], c = 'k', ls = (0,(4,1)), lw = 2, label = 'mean(N=100)')
plt.plot(p_true_range[:20],
RMSER_median[:20], c = 'k', ls = (0,(4,1,2,1)), lw = 2, label =
'median(N=100)')
plt.plot(p_true_range[:20],
RMSER_mode[:20], c = 'k', ls = (0,(1,2)), lw = 2, label = 'MAP(N=100)')
plt.plot(p_true_range[:20],
RMSER_mean2[:20], c = 'k', ls = (0,(4,1)), lw = 1, label = 'mean(N=200)')
plt.plot(p_true_range[:20],
RMSER_median2[:20], c = 'k', ls = (0,(4,1,2,1)), lw = 1, label =
'median(N=200)')
plt.plot(p_true_range[:20],
RMSER_mode2[:20], c = 'k', ls = (2,(2,2)), lw = 1, label = 'MAP(N=200)')
plt.xlabel(r'$\theta_0$', fontsize =
14)
plt.legend(fontsize = 14)
plt.tight_layout()
plt.savefig('Figure4.png')
plt.show()
リスト３ Table 1（1変量正規分布）の計算。リスト４のStanスクリプトを用いている。
import pystan
import numpy as np
import scipy.stats as ss
import matplotlib.pyplot as plt
import seaborn as sb
import pickle
import random
import tkinter as tk
sm = pystan.StanModel(file =
'OneVarNormal.stan')
with open('sm.pkl', 'wb') as f:
pickle.dump(sm, f)
n_simu = 1000
fout = open(f'ResultsCheck.txt', 'w')
Biases_mu_mean = []
Biases_mu_median = []
Biases_mu_map = []
SEBiases_mu_mean = []
SEBiases_mu_median = []
SEBiases_mu_map = []
Biases_sgm_mean = []
Biases_sgm_median = []
Biases_sgm_map = []
SEBiases_sgm_mean = []
SEBiases_sgm_median = []
SEBiases_sgm_map = []
rms_mu_mean = []
rms_mu_median = []
rms_mu_map = []
se_rms_mu_mean = []
se_rms_mu_median = []
se_rms_mu_map = []
rms_sgm_mean = []
rms_sgm_median = []
rms_sgm_map = []
se_rms_sgm_mean = []
se_rms_sgm_median = []
se_rms_sgm_map = []
n_samples = 20
mu_means = []
mu_medians = []
mu_maps = []
sgm_means = []
sgm_medians = []
sgm_maps = []
root = tk.Tk()
root.attributes('-topmost', True)
cnvs = tk.Canvas(root, bg = '#ffffff',
width = 400, height = 200)
cnvs.pack()
cnvs.create_text(200, 100, text =
'Simulation started.',font = ('Gothic', 48), tags = 'mytext')
cnvs.update()
for trial in range(n_simu):
cnvs.itemconfig('mytext',
text = '{0}/{1}'.format(trial+1, n_simu))
cnvs.update()
X = ss.norm().rvs(size =
n_samples)
Data = {'N': len(X), 'X': X}
with open('sm.pkl', 'rb') as
f:
sm =
pickle.load(f)
fit = sm.sampling(data =
Data, control=dict(adapt_delta=0.97, max_treedepth = 30), n_jobs = 1)
Smpl_mu = fit['mu']
Mean_E = Smpl_mu.mean()
mu_means.append(Mean_E)
Med_E =
np.percentile(Smpl_mu, 50)
mu_medians.append(Med_E)
L_mu =
np.percentile(Smpl_mu, 2.5)
U_mu =
np.percentile(Smpl_mu, 97.5)
n_points = 10000
coord = np.linspace(L_mu, U_mu,
n_points)
est_pdf =
ss.gaussian_kde(Smpl_mu).pdf(coord)
map_idx = np.argmax(est_pdf)
MAP_E = coord[map_idx]
mu_maps.append(MAP_E)
Smpl_sgm = fit['sgm']
Mean_E = Smpl_sgm.mean()
sgm_means.append(Mean_E)
Med_E = np.percentile(Smpl_sgm,
50)
sgm_medians.append(Med_E)
L_sgm = 0.0 # Smpl_sgm.min()
U_sgm =
np.percentile(Smpl_sgm, 99)
n_points = 10000
coord = np.linspace(L_sgm,
U_sgm, n_points)
est_pdf =
ss.gaussian_kde(Smpl_sgm).pdf(coord)
map_idx = np.argmax(est_pdf)
MAP_E = coord[map_idx]
sgm_maps.append(MAP_E)
root.destroy()
mu_means = np.array(mu_means)
mu_medians = np.array(mu_medians)
mu_maps = np.array(mu_maps)
mu_mean_bias = mu_means.mean() - 0.0
mu_median_bias = mu_medians.mean() -
0.0
mu_map_bias = mu_maps.mean() - 0
mu_mean_se = ss.sem(mu_means)
mu_median_se = ss.sem(mu_medians)
mu_map_se = ss.sem(mu_maps)
Biases_mu_mean.append(mu_mean_bias)
Biases_mu_median.append(mu_median_bias)
Biases_mu_map.append(mu_map_bias)
SEBiases_mu_mean.append(mu_mean_se)
SEBiases_mu_median.append(mu_median_se)
SEBiases_mu_map.append(mu_map_se)
mss_mu_means = ((mu_means - 0.0) **
2).mean()
rms_mu_mean.append(mss_mu_means **
0.5)
t_sem = ss.sem((mu_means - 0.0) ** 2)
se_rms_mu_mean.append(t_sem ** 0.5)
mss_mu_medians = ((mu_medians - 0.0)
** 2).mean()
rms_mu_median.append(mss_mu_medians **
0.5)
t_sem = ss.sem((mu_medians - 0.0) **
2)
se_rms_mu_median.append(t_sem ** 0.5)
mss_mu_maps = ((mu_maps - 0.0) **
2).mean()
rms_mu_map.append(mss_mu_maps ** 0.5)
t_sem = ss.sem((mu_maps - 0.0) ** 2)
se_rms_mu_map.append(t_sem ** 0.5)
sgm_means = np.array(sgm_means)
sgm_medians = np.array(sgm_medians)
sgm_maps = np.array(sgm_maps)
sgm_mean_bias = sgm_means.mean() - 1.0
sgm_median_bias = sgm_medians.mean() -
1.0
sgm_map_bias = sgm_maps.mean() - 1.0
sgm_mean_se = ss.sem(sgm_means)
sgm_median_se = ss.sem(sgm_medians)
sgm_map_se = ss.sem(sgm_maps)
Biases_sgm_mean.append(sgm_mean_bias)
Biases_sgm_median.append(sgm_median_bias)
Biases_sgm_map.append(sgm_map_bias)
SEBiases_sgm_mean.append(sgm_mean_se)
SEBiases_sgm_median.append(sgm_median_se)
SEBiases_sgm_map.append(sgm_map_se)
mss_sgm_means = ((sgm_means - 0.0) **
2).mean()
rms_sgm_mean.append(mss_sgm_means **
0.5)
t_sem = ss.sem((sgm_means - 0.0) ** 2)
se_rms_sgm_mean.append(t_sem ** 0.5)
mss_sgm_medians = ((sgm_medians - 0.0)
** 2).mean()
rms_sgm_median.append(mss_sgm_medians
** 0.5)
t_sem = ss.sem((sgm_medians - 0.0) **
2)
se_rms_sgm_median.append(t_sem ** 0.5)
mss_sgm_maps = ((sgm_maps - 0.0) **
2).mean()
rms_sgm_map.append(mss_sgm_maps **
0.5)
t_sem = ss.sem((sgm_maps - 0.0) ** 2)
se_rms_sgm_map.append(t_sem ** 0.5)
fout.write('*/n_samples\n')
fout.write(f' {n_samples}')
fout.write('\n')
fout.write('*/Biases_mu_mean\n')
for v in Biases_mu_mean:
fout.write(f' {v}')
fout.write('\n')
fout.write('*/Biases_mu_median\n')
for v in Biases_mu_median:
fout.write(f' {v}')
fout.write('\n')
fout.write('*/Biases_mu_map\n')
for v in Biases_mu_map:
fout.write(f' {v}')
fout.write('\n')
fout.write('*/SEBiases_mu_mean\n')
for v in SEBiases_mu_mean:
fout.write(f' {v}')
fout.write('\n')
fout.write('*/SEBiases_mu_median\n')
for v in SEBiases_mu_median:
fout.write(f' {v}')
fout.write('\n')
fout.write('*/SEBiases_mu_map\n')
for v in SEBiases_mu_map:
fout.write(f' {v}')
fout.write('\n')
fout.write('*/Biases_sgm_mean\n')
for v in Biases_sgm_mean:
fout.write(f' {v}')
fout.write('\n')
fout.write('*/Biases_sgm_median\n')
for v in Biases_sgm_median:
fout.write(f' {v}')
fout.write('\n')
fout.write('*/Biases_sgm_map\n')
for v in Biases_sgm_map:
fout.write(f' {v}')
fout.write('\n')
fout.write('*/SEBiases_sgm_mean\n')
for v in SEBiases_sgm_mean:
fout.write(f' {v}')
fout.write('\n')
fout.write('*/SEBiases_sgm_median\n')
for v in SEBiases_sgm_median:
fout.write(f' {v}')
fout.write('\n')
fout.write('*/SEBiases_sgm_map\n')
for v in SEBiases_sgm_map:
fout.write(f' {v}')
fout.write('\n')
fout.write('*/rms_mu_mean\n')
for v in rms_mu_mean:
fout.write(f' {v}')
fout.write('\n')
fout.write('*/rms_mu_median\n')
for v in rms_mu_median:
fout.write(f' {v}')
fout.write('\n')
fout.write('*/rms_mu_map\n')
for v in rms_mu_map:
fout.write(f' {v}')
fout.write('\n')
fout.write('*/se_rms_mu_mean\n')
for v in se_rms_mu_mean:
fout.write(f' {v}')
fout.write('\n')
fout.write('*/se_rms_mu_median\n')
for v in se_rms_mu_median:
fout.write(f' {v}')
fout.write('\n')
fout.write('*/se_rms_mu_map\n')
for v in se_rms_mu_map:
fout.write(f' {v}')
fout.write('\n')
fout.write('*/rms_sgm_mean\n')
for v in rms_sgm_mean:
fout.write(f' {v}')
fout.write('\n')
fout.write('*/rms_sgm_median\n')
for v in rms_sgm_median:
fout.write(f' {v}')
fout.write('\n')
fout.write('*/rms_sgm_map\n')
for v in rms_sgm_map:
fout.write(f' {v}')
fout.write('\n')
fout.write('(*/se_rms_sgm_mean\n')
for v in se_rms_sgm_mean:
fout.write(f' {v}')
fout.write('\n')
fout.write('*/se_rms_sgm_median\n')
for v in se_rms_sgm_median:
fout.write(f' {v}')
fout.write('\n')
fout.write('*/se_rms_sgm_map\n')
for v in se_rms_sgm_map:
fout.write(f' {v}')
fout.write('\n')
fout.close()
data {
int N;
real X[N];
}
parameters {
real mu;
real<lower = 0.0> sgm;
}
model {
mu ~ normal(0.0, 1000.0);
sgm ~ exponential(0.001);
for (i in 1:N) {
X[i]
~ normal(mu, sgm);
}
}
リスト５ Tables 2, 3（相関係数）の計算。リスト６のStanスクリプトを用いている。
import pystan
import numpy as np
import scipy.stats as ss
import matplotlib.pyplot as plt
import seaborn as sb
import pickle
import random
import tkinter as tk
sm = pystan.StanModel(file =
'TwoVarNormal.stan')
with open('sm.pkl', 'wb') as f:
pickle.dump(sm, f)
n_samples = 50
n_simu = 200
rhos = [-0.9, -0.6, -0.3, 0.0, 0.3,
0.6, 0.9]
bias_mean = []
bias_median = []
bias_MAP = []
SE_bias_mean = []
SE_bias_median = []
SE_bias_MAP = []
RMSE_mean = []
RMSE_median = []
RMSE_MAP = []
SE_RMSE_mean = []
SE_RMSE_median = []
SE_RMSE_MAP = []
root = tk.Tk()
root.attributes('-topmost', True)
cnvs = tk.Canvas(root, bg = '#ffffff',
width = 400, height = 300)
cnvs.pack()
cnvs.create_text(200, 150, text =
'Simulation started.',font = ('Gothic', 32), tags = 'mytext')
cnvs.update()
for rho_true in rhos:
#print('rho =',
rho_true)
s_bias_mean = 0.0
s_bias_median = 0.0
s_bias_MAP = 0.0
se_s_bias_mean = []
se_s_bias_median = []
se_s_bias_MAP = []
ss_mean = 0.0
ss_median = 0.0
ss_MAP = 0.0
se_ss_mean = []
se_ss_median = []
se_ss_MAP = []
Means = []
Medians = []
MAPs = []
for trial in range(n_simu):
cnvs.itemconfig('mytext', text = f'rho = {rho_true}\ntrial =
{trial+1}/{n_simu}')
cnvs.update()
X =
0
u =
np.random.normal(size = n_samples)
v =
np.random.normal(size = n_samples)
random.shuffle(u)
random.shuffle(v)
t =
u * ((1.0 - (rho_true**2))**0.5) + v * rho_true
X =
np.vstack([t, v]).T
Data
= {'N': len(X), 'X': X}
sm =
0
with
open('sm.pkl', 'rb') as f:
sm = pickle.load(f)
fit
= 0
fit
= sm.sampling(data = Data, iter =
4000, n_jobs = 1)
Smpl_rho = fit['rho']
Mean_E = Smpl_rho.mean()
Med_E = np.percentile(Smpl_rho, 50)
coord = np.linspace(-1.0, 1.0, 20000)
est_pdf = ss.gaussian_kde(Smpl_rho).pdf(coord)
map_idx = np.argmax(est_pdf)
MAP_E = coord[map_idx]
Means.append(Mean_E)
Medians.append(Med_E)
MAPs.append(MAP_E)
s_bias_mean += Mean_E - rho_true
s_bias_median += Med_E - rho_true
s_bias_MAP += MAP_E - rho_true
se_s_bias_mean.append(Mean_E - rho_true)
se_s_bias_median.append(Med_E
- rho_true)
se_s_bias_MAP.append(MAP_E - rho_true)
ss_mean += (Mean_E - rho_true) ** 2
ss_median += (Med_E - rho_true) ** 2
ss_MAP += (MAP_E - rho_true) ** 2
se_ss_mean.append((Mean_E - rho_true) ** 2)
se_ss_median.append((Med_E - rho_true) ** 2)
se_ss_MAP.append((MAP_E - rho_true) ** 2)
bias_mean.append(s_bias_mean
/ n_simu)
bias_median.append(s_bias_median / n_simu)
bias_MAP.append(s_bias_MAP /
n_simu)
SE_bias_mean.append(ss.sem(se_s_bias_mean))
SE_bias_median.append(ss.sem(se_s_bias_median))
SE_bias_MAP.append(ss.sem(se_s_bias_MAP))
RMSE_mean.append((ss_mean /
n_simu) ** 0.5)
RMSE_median.append((ss_median / n_simu) ** 0.5)
RMSE_MAP.append((ss_MAP / n_simu) ** 0.5)
SE_RMSE_mean.append(ss.sem(se_ss_mean) ** 0.5)
SE_RMSE_median.append(ss.sem(se_ss_median) ** 0.5)
SE_RMSE_MAP.append(ss.sem(se_ss_MAP) ** 0.5)
root.destroy()
fout =
open(f'ResultsN{n_samples}.txt', 'w')
print()
print('bias_mean =', bias_mean)
print('bias_edian =', bias_median)
print('bias_MAP =', bias_MAP)
print('RMSE_mean =', RMSE_mean)
print('RMSE_median =', RMSE_median)
print('RMSE_MAP =', RMSE_MAP)
fout.write('*/bias_mean\n')
for v in bias_mean:
fout.write(f' {v}')
fout.write('\n')
fout.write('*/SE_bias_mean\n')
for v in SE_bias_mean:
fout.write(f' {v}')
fout.write('\n')
fout.write('*/bias_median\n')
for v in bias_median:
fout.write(f' {v}')
fout.write('\n')
fout.write('*/SE_bias_median\n')
for v in SE_bias_median:
fout.write(f' {v}')
fout.write('\n')
fout.write('*/bias_MAP\n')
for v in bias_MAP:
fout.write(f' {v}')
fout.write('\n')
fout.write('*/SE_bias_MAP\n')
for v in SE_bias_MAP:
fout.write(f' {v}')
fout.write('\n')
fout.write('*/RMSE_mean\n')
for v in RMSE_mean:
fout.write(f' {v}')
fout.write('\n')
fout.write('*/SE_RMSE_mean\n')
for v in SE_RMSE_mean:
fout.write(f' {v}')
fout.write('\n')
fout.write('*/RMSE_median\n')
for v in RMSE_median:
fout.write(f' {v}')
fout.write('\n')
fout.write('*/SE_RMSE_median\n')
for v in SE_RMSE_median:
fout.write(f' {v}')
fout.write('\n')
fout.write('*/RMSE_MAP\n')
for v in RMSE_MAP:
fout.write(f' {v}')
fout.write('\n')
fout.write('*/SE_RMSE_MAP\n')
for v in SE_RMSE_MAP:
fout.write(f' {v}')
fout.write('\n')
fout.close()
data {
int N;
vector[2] X[N];
}
parameters {
real<lower = -1.0, upper
= 1.0> rho;
real mu1;
real mu2;
real<lower = 0.0>
sgm1;
real<lower = 0.0> sgm2;
}
transformed parameters {
vector[2] mu;
cov_matrix[2] cov;
mu[1] = mu1;
mu[2] = mu2;
cov[1][1] = sgm1 * sgm1;
cov[1][2] = rho * sgm1 *
sgm2;
cov[2][1] = rho * sgm1 *
sgm2;
cov[2][2] = sgm2 * sgm2;
}
model {
rho ~ uniform(-1.0, 1.0);
mu1 ~ normal(0.0, 1000.0);
mu2 ~ normal(0.0, 1000.0);
sgm1 ~ exponential(0.001);
sgm2 ~ exponential(0.001);
for (i in 1:N) {
X[i]
~ multi_normal(mu, cov);
}
}
MAP推定値の算出スクリプト例
ベイズ分析（岡本、2022、pp.160-167）では、分析結果は事後分布として得られる。パラメータの点推定値は、事後分布からMAP推定値（岡本、2019、pp.175-176）として求めることができる。例えば、次の関数calcMAPEstのように事後分布のサンプルsamplesから算出できる。
def calcMAPEst(samples, a = 0.05,
n_points = 1000):
""" MAP推定値とMAP推定値における推定確率密度関数の値を求める
samples: MAP推定値を求めるサンプル
a: サンプルの左右a/2を除く。デフォルト値=0.05
n_points: KDEによる確率密度関数値を求める点の数。
デフォルト値=1000
"""
import numpy as np
import scipy.stats as ss
Lp, Up = np.percentile(samples,
[100 * a/2, 100 * (1 - a/2)])
coord = np.linspace(Lp, Up,
n_points)
est_pdf =
ss.gaussian_kde(samples).pdf(coord)
map_idx = np.argmax(est_pdf)
MAP_Est =
coord[map_idx]
return MAP_Est,
est_pdf[map_idx] # MAP推定値とMAP推定値における確率密度関数値
上のスクリプトでは、サンプルの両端の部分を切り捨てている。これは、サンプルの裾が以上に長い場合に、モード付近の確率密度分布関数の精度が落ちないようにするためである。サンプルから両端を切り取った残りの区間をn_points個の点で離散化して、KDE (Kernel Density Estimation) により確率密度分布関数の推定を行っている。これらn_points個の点における値の最大値を与える点として、MAP推定値を求め、その点における確率密度関数の値とともに関数calcMAPEstの返り値としている。
MAP推定値を求める例として、リストA1に正規分布のパラメータ平均値と標準偏差を求めるスクリプト例を示す。
リスト１A 正規分布のパラメータ平均値と標準偏差を求めるスクリプト例
import numpy as np
import matplotlib.pyplot as plt
import seaborn as sb
import pystan
import scipy.stats as ss
def calcMAPEst(samples, a = 0.05,
n_points = 1000):
""" MAP推定値とMAP推定値における推定確率密度関数の値を求める
samples: MAP推定値を求めるサンプル
a: サンプルの左右a/2を除く。デフォルト値=0.05
n_points: KDEによる確率密度関数値を求める点の数。
デフォルト値=1000
"""
import numpy as np
import scipy.stats as ss
Lp, Up =
np.percentile(samples, [100 * a/2, 100 * (1 - a/2)])
coord = np.linspace(Lp, Up,
n_points)
est_pdf =
ss.gaussian_kde(samples).pdf(coord)
map_idx = np.argmax(est_pdf)
MAP_Est =
coord[map_idx]
return MAP_Est,
est_pdf[map_idx] # MAP推定値とMAP推定値における確率密度関数値
Data = [51, 38, 60, 42, 47]
#
#
正規分布モデル
#
StanCode_normal = """
data {
int
n;
real
d[n];
}
parameters {
real
mu;
real<lower = 0.0> sigma;
}
model {
mu ~
normal(0.0, 1000.0);
sigma ~ exponential(0.001);
d ~
normal(mu, sigma);
}
"""
sm = pystan.StanModel(model_code =
StanCode_normal)
fit = sm.sampling(data =
{'n':len(Data), 'd':Data}, iter = 10000, n_jobs = 1)
print(fit)
mu_MAP = calcMAPEst(fit['mu'])[0]
sigma_MAP =
calcMAPEst(fit['sigma'])[0]
#
# グラフ描画
#
plt.figure(figsize = (8,3))
#
# muの事後分布の描画
#
plt.subplot(121)
sb.kdeplot(fit['mu'])
plt.xlabel('$\mu$', fontsize = 14)
plt.title('$\hat\mu =$
{0:.1f}(MAP)'.format(mu_MAP), fontsize = 14)
#
#
sigmaの事後分布の描画
#
plt.subplot(122)
sb.kdeplot(fit['sigma'])
plt.xlabel('$\sigma$', fontsize = 14)
plt.title('$\hat\sigma =$
{0:.1f}(MAP)'.format(sigma_MAP), fontsize = 14)
plt.tight_layout()
plt.savefig('PostMuSgm.png') # グラフの保存
plt.show()
print('PostMuSgm.png was saved.')
リストA1のスクリプトを実行すると、Stan（岡本、2019、pp.179-210）によるMCMCサンプリングが行われる。その結果は、図A1のように表示される。
図A1
MAP推定値がグラフの上に表示されている。
図A1の表示されているウィンドウを閉じると、スクリプトの実行終了となる。端末には、以下のメッセージが表示されている。
Inference for Stan model:
anon_model_b3c1b2700ba7728d2ab9aa18f53eeba1.
4 chains, each with iter=10000;
warmup=5000; thin=1;
post-warmup draws per chain=5000,
total post-warmup draws=20000.
mean
se_mean sd 2.5% 25% 50% 75% 97.5% n_eff Rhat
mu 47.61 0.15 7.76 33.6 44.33 47.61 50.89 61.07 2565 1.0
sigma 13.34 0.18 9.53 5.51 8.28 10.91 15.25 35.84 2785 1.0
lp__ -11.86 0.02 1.31 -15.38 -12.36 -11.46
-10.93 -10.6 3396 1.0
Samples were drawn using NUTS at Thu
Dec 15 08:36:13 2022.
For each parameter, n_eff is a crude
measure of effective sample size,
and Rhat is the potential scale
reduction factor on split chains (at
convergence, Rhat=1).
PostMuSgm.png was saved.
平均値のパラメータmuの推定値は、事後分布の平均値、中央値、MAP推定値に大差はないが、これは事後分布がほぼ左右対称形であることによる。事後分布が左右対称形でないsigmaの場合は、平均値13.34、中央値10.91、MAP推定値8.8と異なった値になっている。
なお、実行は、PyStan2.17.1.0/Python3.7/Anaconda/Windows 10で行った。
PyStanについては、＜岡本安晴「いまさら聞けないPythonでデータ分析」丸善出版＞pp.178-210で説明している。