拙著第９章「確率計算とPythonスクリプト」補足説明
MAP推定値の勧め
モデルのパラメータの点推定値をベイズ分析で求めるときは、事後分布から推定することになる。分布の代表値として、平均値、中央値、モードがあるが、どれを点推定値として採用するのがよいかの問題がある。Gelmanら（2021）では、分布の代表値として平均値あるいは中央値が用いられるが、中央値の方が計算が安定しているとして、中央値に基づく方法が採用されている（p. 73）。モードについては、Gelmanら（2021）は、データと事前分布に対して最良の適合を与えるものであると説明している（p. 107）。
パラメータの真値が分かっているシミュレーションにおいて、上下法データの場合、３つの推定値の中でモード（MAP推定値）が最も良かったことが報告されている（岡本、2011、Table 1）。多件法による上下法の場合のシミュレーションの場合も、MAP推定値が中央値よりよかったと報告されている（岡本、2019、Tables 2,3,4）。上下法以外の場合を調べたものとして、比率、平均値と標準偏差、相関係数を扱ったものがあるが、3つの推定値の中でMAP推定値が良かったと報告されている（岡本、2022）。
３つの点推定値の比較例を、平均値と標準偏差についての簡単なシミュレーション例で示す。
使用したスクリプトファイルは、圧縮ファイルc09files.zipにまとめた。ダウンロード展開して自由に利用されたい。
正規分布モデルによってベイズ分析したときの事後分布は、下図のようになる。
図１a
図１aは、平均値パラメータ
の事後分布である。この場合、３つの代表値、平均値、中央値、モード、はほぼ一致している。
しかし、下図１bの場合の分布は対称形ではない。
図１b
図１bは、標準偏差パラメータ
の事後分布である。パラメータの場合は、点推定値として平均値をとるか、中央値をとるか、モードをとるかによって、値は実質的に異なる。
図１aと図１bは、リスト１aとリスト１bのスクリプトcheck_int_s1.pyとcheck_ind_s2.pyを順番に実行すると得ることができる。リスト１aのスクリプトはpystan 3を利用しているが、Windowsで実行する場合は、WSL（Windows Subsystem for Linux）のUbuntuで実行する。Ubuntu/WSL/WindowsにインストールされたAnacondaにPython 3.11仮想環境を構築してPyStan 3をインストールすればよい。詳しくは、このウェブサイトを参照されたい。
３つの点推定値の比較のためのスクリプトをリスト２a（check_est_s1.py）およびリスト２b（check_est_s2.py）に用意した。リスト２aのスクリプトを実行してから、リスト２bのスクリプトを実行すると、下図２aおよび2ｂを得る。この場合も、Windowsの場合は、リスト２aをAnaconda/Ubuntu/WSL/WindowsのPython 3.11仮想環境で実行した後、Windowsの環境でリスト２bを実行する。
図２a
図２aは、平均値が
、標準偏差が
の正規分布からサンプル数
のデータをとりベイズ分析を行うという試行をシミュレーションで100試行行ったときのパラメータの３つの推定値
の分布をKDE密度関数で表したものである。図２aに示されているの３つの推定値の分布はほとんど重なっており、RMS誤差も違いがないと言える。これは、の事後分布が、図１aに示されているように対称形であることによる。
図２b
図２bは、標準偏差の点推定値
の３つのもの、平均値（mean）、中央値（median）、モード（map）、の分布であり、違いが見られる。RMSエラーも平均値、中央値、MAP推定値の順に小さくなっている。
Gelman, A., Hill, J., & Vehtari, A. (2021). Regression and other stories. Cambridge University Press.
岡本安晴（2011）上下法データのベイズ的分析 : 確率モデルによるデータの活用．基礎心理学研究、30、44-55．DOI https://doi.org/10.14947/psychono.KJ00007667868
岡本安晴（2019）多件法による上下法．基礎心理学研究、38、90-104．DOI https://doi.org/10.14947/psychono.38.2
岡本安晴（2022）ベイズ分析における点推定値例のRMSEとbiasによる比較――実験心理学の観点からの検討――．基礎心理学研究、41、1-7．DOI https://doi.org/10.14947/psychono.41.1
リスト１a．１変量正規分布のサンプル例スクリプト―１（check_ind_s1.py）
import numpy
as np
import scipy.stats
as ss
import stan
import pickle
import arviz
as az
import pandas as pd
pd.options.display.max_rows = 1000
pd.options.display.max_columns = 1000
stan_script = """
data {
int N;
array[N] real X;
}
parameters {
real
mu;
real<lower=0.0> sgm;
}
model {
mu ~
normal(.0, 10.0);
sgm ~ exponential(0.1);
X ~
normal(mu, sgm);
}
"""
N = 10
X = ss.norm.rvs(size
= N)
sm = stan.build(stan_script,data = {'N':len(X),
'X':X} )
fit = sm.sample(num_samples=10000)
inference_data = az.from_pystan(posterior
= fit, posterior_model = sm)
smry = az.summary(inference_data)
print('Summary =\n',smry)
df = fit.to_frame()
with open('d_frame.pkl',
'wb') as f:
pickle.dump(df, f)
リスト１b．１変量正規分布のサンプル例スクリプト―２（check_ind_s2.py）
import numpy
as np
import scipy.stats
as ss
import matplotlib.pyplot
as plt
import seaborn
as sb
import pickle
import arviz
as az
import pandas as pd
with open('d_frame.pkl',
'rb') as f:
df
= pickle.load(f)
sb.kdeplot(df['mu'])
plt.xlabel('$\mu$', fontsize
= 16)
plt.savefig('FigMu.png')
plt.show()
sb.kdeplot(df['sgm'])
plt.xlabel('$\sigma$', fontsize
= 16)
plt.savefig('FigSgm.png')
plt.show()
リスト２a．正規分布の場合のシミュレーションスクリプト例―１（check_est_s1.py）
import numpy
as np
import scipy.stats
as ss
import stan
import pickle
import arviz
as az
import pandas as pd
pd.options.display.max_rows = 1000
pd.options.display.max_columns = 1000
stan_script = """
data {
int N;
array[N] real X;
}
parameters {
real
mu;
real<lower=0.0> sgm;
}
model {
mu ~
normal(.0, 10.0);
sgm ~ exponential(0.1);
X ~
normal(mu, sgm);
}
"""
def CalcMAPEst(samples,
a = 0.05, n_points = 2000):
"""
Calculatte a MAP estimate from a KDE graph on [Lp, Up]
Lp and Up are 100*a/2 and 100(1-a/2) percentile points of
samples
"""
Lp,
Up = np.percentile(samples, [100 * a/2, 100 * (1 -
a/2)]) # import numpy
as np
coord
= np.linspace(Lp, Up, n_points)
est_pdf
= ss.gaussian_kde(samples).pdf(coord) # import scipy.stats
as ss
map_idx
= np.argmax(est_pdf)
MAP_Est
= coord[map_idx]
return MAP_Est,
est_pdf[map_idx]
N = 10
n_simu = 100
mu_means = []
mu_medians = []
mu_maps = []
sgm_means = []
sgm_medians = []
sgm_maps = []
ws = ''
for i in
range(90):
ws
+= ' '
for simu in
range(n_simu):
print(ws, f'simu = {simu}/{n_simu} ', end =
'\r')
X = ss.norm.rvs(size
= N)
sm
= stan.build(stan_script,data
= {'N':len(X), 'X':X} )
fit = sm.sample(num_samples=10000)
mu_mean
= np.mean(fit['mu'])
mu_median
= np.median(fit['mu'])
mu_map
= CalcMAPEst(fit['mu'])[0]
mu_means.append(mu_mean)
mu_medians.append(mu_median)
mu_maps.append(mu_map)
sgm_mean
= np.mean(fit['sgm'])
sgm_median
= np.median(fit['sgm'])
sgm_map
= CalcMAPEst(fit['sgm'])[0]
sgm_means.append(sgm_mean)
sgm_medians.append(sgm_median)
sgm_maps.append(sgm_map)
with open('temp.pkl',
'wb') as f:
pickle.dump((mu_means, mu_medians, mu_maps,
sgm_means, sgm_medians,
sgm_maps), f)
print('temp.pkl
was saved.')
リスト２b．正規分布の場合のシミュレーションスクリプト例―２（check_est_s2.py）
import numpy
as np
import scipy.stats
as ss
import matplotlib.pyplot
as plt
import seaborn
as sb
import pickle
import arviz
as az
import pandas as pd
def CalcMAPEst(samples,
a = 0.05, n_points = 2000):
"""
Calculatte a MAP estimate from a KDE graph on [Lp, Up]
Lp and Up are 100*a/2 and 100(1-a/2) percentile points of
samples
"""
Lp,
Up = np.percentile(samples, [100 * a/2, 100 * (1 -
a/2)]) # import numpy
as np
coord
= np.linspace(Lp, Up, n_points)
est_pdf
= ss.gaussian_kde(samples).pdf(coord) # import scipy.stats
as ss
map_idx
= np.argmax(est_pdf)
MAP_Est
= coord[map_idx]
return MAP_Est,
est_pdf[map_idx]
with open('temp.pkl',
'rb') as f:
(mu_means,
mu_medians, mu_maps, sgm_means, sgm_medians, sgm_maps) = \
pickle.load(f)
mu_means = np.array(mu_means)
mu_medians = np.array(mu_medians)
mu_maps = np.array(mu_maps)
rms_mu_means = np.mean((mu_means - 0)**2)
rms_mu_medians = np.mean((mu_medians - 0)**2)
rms_mu_maps = np.mean((mu_maps - 0)**2)
print('RMS(mu_mean)
=', rms_mu_means)
print('RMS(mu_median)
=', rms_mu_medians)
print('RMS(mu_map)
=', rms_mu_maps)
sgm_means = np.array(sgm_means)
sgm_medians = np.array(sgm_medians)
sgm_maps = np.array(sgm_maps)
rms_sgm_means = np.mean((sgm_means - 1)**2)
rms_sgm_medians = np.mean((sgm_medians - 1)**2)
rms_sgm_maps = np.mean((sgm_maps - 1)**2)
print('\nRMS(sgm_mean) =', rms_sgm_means)
print('RMS(sgm_median)
=', rms_sgm_medians)
print('RMS(sgm_map)
=', rms_sgm_maps)
sb.kdeplot(mu_means,
label = 'mean')
sb.kdeplot(mu_medians,
label = 'median')
sb.kdeplot(mu_maps,
label = 'map')
plt.plot([0, 0], [0, 1.2], c = 'r', ls = '--',
label ='$\mu = 0$')
plt.title(f'RMSEs of
estimates of $\mu$\n{rms_mu_means:.3f}(for means)' +
f', {rms_mu_medians:.3f}(for medians))' +
f', {rms_mu_maps:.3f}(for
maps)', fontsize = 14)
plt.xlabel('$\hat\mu$', fontsize
= 16)
plt.legend()
plt.tight_layout()
plt.savefig('FigMuHat.png')
plt.show()
sb.kdeplot(sgm_means,
label = 'mean')
sb.kdeplot(sgm_medians,
label = 'median')
sb.kdeplot(sgm_maps,
label = 'map')
plt.plot([1, 1], [0, 1.4], c = 'r', ls = '--',
label = '$\sigma = 1$')
plt.title(f'RMSEs of
estimates of $\sigma$\n{rms_sgm_means:.3f}(for means)' +
f', {rms_sgm_medians:.3f}(for
medians))' +
f', {rms_sgm_maps:.3f}(for
maps)', fontsize = 14)
plt.xlabel('$\hat\sigma$', fontsize
= 16)
plt.legend()
plt.savefig('FigSgmHat.png')
plt.show()