ベイズ点推定値と最高密度区間
事後分布に対する代表値として平均値、中央値、mode（MAP推定値）が挙げられる。これらの点推定値を、最高密度区間（Highest Density Interval: HDI）と比較してみた。事後分布の推定法としてMCMCサンプリングがあり、平均値と中央値はサンプルの平均値と中央値によって与えられる。サンプルは離散値であるので、事後分布が連続量のとき、modeをサンプル値のヒストグラムから最頻値として算出することには無理がある。ヒストグラムからではなく、事後分布の連続量としての分布をサンプルからKDE（kernel density estimation）により推定して、mode（MAP推定値）を求める求めることができる。このPythonスクリプトを以下のように作成した。
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]
上のスクリプトでは、MCMCによるサンプルsamplesからMAP推定値を求めMAP_Estで表している。
最奥密度区間HDIはChen-Shaoのアルゴリズムを利用して求めた。
点推定値とHDIの比較のため、次のデータ
data = [57,
30, 34, 59, 63, 26, 60, 58, 34, 56]
に対して、正規分布モデルを適用してパラメータ平均値
と標準偏差
の事後分布を求めるStanスクリプトを利用するPythonスクリプトを以下のように作成した。
UniVarNormStan = """
data {
int
N;
real Data[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) {
Data[i] ~ normal(mu, sgm);
}
}
"""
data = [57, 30, 34, 59, 63, 26, 60,
58, 34, 56]
sm = pystan.StanModel(model_code =
UniVarNormStan)
fit = sm.sampling(data = {'Data':data,
'N':len(data)}, iter = 4000, n_jobs
= 1)
上のスクリプトを含むリスト１のプログラムを実行すると図１の結果を得た。
図１
平均値の事後分布はほぼ左右対称形であるので、３つの代表値はほぼ一致している。しかし、標準偏差の事後分布は裾が正の方向に伸びているので、３つの代表値は異なり、MAP推定値、中央値、平均値の順に大きくなっている。
岡本（2011）は、３つの点推定値の比較を心理物理学実験のシミュレーションにおいて行っている。シミュレーションであるのでパラメータの真値は分かっており、真値の点推定値の比較をbiasとRMS errorで行っている。PSEについては、事後分布は左右対称形であり、３つの点推定値に差は認められなかった。JNDについては、分布は非対称形で、MAP推定値が最も良い推定値であり、平均が最も悪い推定値であったと報告されている（Table 1）。
上図１の下段に、50％最高密度区間（50％HDI）と３つの推定値との関係が表されている。MAP推定値はいずれもHDIの中ほどにあるが、平均値の方は、σの事後分布の場合、50％HDIの右端に位置していて、HDIから外れそうである。
Gelman, A., Hill, J., & Vehtari, A. (2021).”Regression and Other Stories.” Cambridge Univ. Pressのp. 73には
We typically prefer median-based summaries because they are more computationally stable
とあり、p. 107には
For a Bayesian model, the corresponding point estimate is the posterior mode, which provides the best overall fit to data and prior distribution.
とある。ここで、modeとはMAP推定値であることに注意。
本ウエブサイトで使用したプログラムの完全なリストをリスト１に示す。リスト１のプログラムは圧縮ファイルcalcmapestfile.zipにまとめた。
リスト１ MAP推定値算出スクリプト例
import numpy as np
import pystan
import matplotlib.pyplot as plt
import seaborn as sb
import scipy.stats as ss
UniVarNormStan = """
data {
int
N;
real
Data[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) {
Data[i] ~ normal(mu, sgm);
}
}
"""
data = [57, 30, 34, 59, 63, 26, 60,
58, 34, 56]
sm = pystan.StanModel(model_code =
UniVarNormStan)
fit = sm.sampling(data = {'Data':data,
'N':len(data)},
n_jobs = 1)
# The parameter n_jobs
= 1 is for Windows
print(fit)
def CalcMAPEst(samples, a = 0.05,
n_points = 10000):
"""
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]
def calcHDI(sample, a):
"""
Calculation of 100(1-a)% HDI for sample
Chen-Shao HPD Estimation algorithm
"""
smpl = np.sort(sample)
n = len(smpl)
ckW = max(smpl) - min(smpl)
cki = 0
RB = int(n * (1-a))
L0 = 0
U0 = n-1
while (RB + cki < n):
if
smpl[RB + cki] - smpl[cki] < ckW:
ckW =
smpl[RB + cki] - smpl[cki]
L0 = cki
U0 = RB + cki
cki
+= 1
print('L0 =', L0, ' U0 =', U0)
ckW = max(smpl) - min(smpl)
cki = 1
LB = n - int(n * (1-a))
L1 = 0
U1 = n-1
while (LB - cki >= 0):
if
smpl[n - cki] - smpl[LB - cki] < ckW:
ckW = smpl[n - cki] - smpl[LB - cki]
U1 = n - cki
L1 = LB - cki
cki
+= 1
print('L1 =', L1, ' U1 =', U1)
return smpl[int((L1+L0)/2)],
smpl[int((U1+U0)/2)]
fig, ax = plt.subplots(figsize = (12,
6))
mean_mu = fit['mu'].mean()
med_mu = np.percentile(fit['mu'], 50)
map_mu, max_v = CalcMAPEst(fit['mu'])
print(f'MAP_E for mu = {map_mu:.2f}')
L, U = calcHDI(fit['mu'], 0.5)
print(f'50% HDI = [{L:.3f}, {U:.3f}]')
plt.subplot(2,2,1)
plt.title('Posterior Distribution of
$\mu$')
sb.kdeplot(fit['mu'],
clip=np.percentile(fit['mu'], [2.5, 97.5]))
plt.plot([map_mu, map_mu], [0.0,
max_v], label = f'MAP_Est({map_mu:.2f})')
plt.plot([med_mu, med_mu], [0.0,
max_v], label = f'Median({med_mu:.2f})')
plt.plot([mean_mu, mean_mu], [0.0,
max_v], label = f'Mean({mean_mu:.2f})')
plt.legend()
plt.subplot(2,2,3)
plt.title('50% HDI, Mean, Median and
MAP est.\n' + 'of Posterior Distribution of $\mu$')
plt.plot([L, U], [0.0, 0.0], lw = 5,
color = 'y', label = '50% HDI')
plt.plot(map_mu, 0.0, marker = 's',
markersize = 10, alpha = 0.5, label = 'MAP')
plt.plot(med_mu, 0.0, marker = '^',
markersize = 10, alpha = 0.5, label = 'Median')
plt.plot(mean_mu, 0.0, marker = 'v',
markersize = 10, alpha = 0.5, label = 'Mean')
plt.yticks([])
plt.tight_layout()
plt.legend()
mean_sgm = fit['sgm'].mean()
med_sgm = np.percentile(fit['sgm'],
50)
map_sgm, max_v =
CalcMAPEst(fit['sgm'])
print(f'MAP_E for sigma = {map_sgm:.2f}')
L, U = calcHDI(fit['sgm'], 0.5)
print(f'50% HDI = [{L:.3f}, {U:.3f}]')
plt.subplot(2,2,2)
plt.title('Posterior Distribution of
$\sigma$')
sb.kdeplot(fit['sgm'],
clip=np.percentile(fit['sgm'], [2.5, 97.5]))
plt.plot([map_sgm, map_sgm], [0.0, max_v],
label = f'MAP_Est({map_sgm:.2f})')
plt.plot([med_sgm, med_sgm], [0.0,
max_v], label = f'Median({med_sgm:.2f})')
plt.plot([mean_sgm, mean_sgm], [0.0,
max_v], label = f'Mean({mean_sgm:.2f})')
v_kde =
ss.gaussian_kde(fit['sgm']).evaluate(np.linspace(L, U, 10000))
plt.legend()
plt.subplot(2,2,4)
plt.title('50% HDI, Mean, Median and
MAP est.\n' + 'of Posterior Distribution of $\sigma$')
plt.plot([L, U], [0.0, 0.0], lw = 5,
color = 'y', label = '50% HDI')
plt.plot(map_sgm, 0.0, marker = 's',
markersize = 10, alpha = 0.5, label = 'MAP')
plt.plot(med_sgm, 0.0, marker = '^',
markersize = 10, alpha = 0.5, label = 'Median')
plt.plot(mean_sgm, 0.0, marker = 'v',
markersize = 10, alpha = 0.5, label = 'Mean')
plt.yticks([])
plt.legend()
plt.tight_layout()
plt.savefig('Figure.png')
plt.show()
なお、PyStan（Python用Stan）の説明を、＜岡本安晴「いまさら聞けないPythonでデータ分析」丸善出版＞で行っている。