Bayesian Point Estimates and Highest Density Interval
Yasuharu Okamoto, 2021.07, 2021.08
As point estimates for a posterior distribution, three estimates are well-known, mean, median, and mode. The mode in this case is called MAP (maximum a posteriori) estimate. As an interval estimate of a parameter, highest density interval (HDI) is well-known.
In case of MCMC sampling, mean and median can be estimated by the sample mean and median. But, sample is consisted of discrete values, the mode cannot be calculated directly from the sample. KDE (kernel density estimate) gives us approximation to the continuous posterior distribution. The following code calculate MAP estimate from a sample.
Lp, Up = np.percentile(sample, [100 * a/2, 100 *
(1 - a/2)]) # import numpy as np
coord = np.linspace(Lp, Up,
n_points)
est_pdf =
ss.gaussian_kde(sample).pdf(coord)
# import scipy.stats as ss
map_idx =
np.argmax(est_pdf)
MAP_Est =
coord[map_idx]
Highest density interval is calculated by Chan-Shao HPD estimation algorithm in Listing 1 below.
To demonstrate relations between HDI and the three estimates, the following data was prepared.
data = [57,
30, 34, 59, 63, 26, 60, 58, 34, 56]
To analyze the data by univariate normal model, the following Stan script was constructed.
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);
}
}
"""
Python script, which analyze the data by the above Stan script is as follows.
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)
# The parameter n_jobs
= 1 is for Windows
Complete Python script is shown in Listing 1. The script file in Listing 1 is included in the archived file calcmaphdifiles.zip , which can be freely used.
Run the script in Listing 1, we get Figure 1.
Figure 1 Mean, median, MAP estimate, and 50% HDI.
For the parameter mu, the posterior distribution is symmetric, so the three estimates are essentially the same. However, the posterior distribution of the parameter sigma is positively skewed, and the three estimates are clearly distinct. In the bottom row, 50% HDIs are shown. The MAP estimates are about at the center of the HDIs. However, the mean of the posterior distribution of parameter sigma is near the right edge of the HDI. Which of the three estimates is the best one?
Lambert (2018) says that the mean is the best. Gelman et al. (2021) recommend the median, because of stability. These discussions are on representative values of a distribution.
However, when we consider relations between the true values and estimates, we get another view. Table 1 in Okamoto (2011) shows that the MAP estimate is the best, and the mean is the worst. Okamoto (2011) conducted psychophysical simulation, in which PSEs and JNDs were estimated. The parameters corresponding to PSE and JND were denoted by mu and sigma. In simulation, true values of the parameters are known. For the parameter mu (PSE), three estimates, mean, median, and MAP estimate, were essentially the same with respect to biases and RMS errors, but for the parameter sigma (JND), MAP estimates were the best and means the worst.
Lambert (2018) discusses that MAP estimates should not be used, because for a bimodal distribution, the mode (MAP estimate) may not reflect the main body of the distribution. However, in this bimodal case, two MAP estimates, each of which represents the respective mode, should be calculated, instead of averaging the whole distribution. Each mode may suggests the respective component of the distribution.
Gelman, A., Hill, J., & Vehtari, A. (2021). Regression and other stories. Cambridge University Press.
Lambert, B. (2018). A student’s guide to Bayesian statistics. SAGE.
Okamoto, Y. (2011). Bayesian analysis of data by the up-down method : Efficient use of data with a stochastic model. The Japanese Journal of Psychonomic Science, 30, 44-55. https://doi.org/10.14947/psychono.KJ00007667868
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()