Estimation of the mean, median, and mode of a
posterior distribution from MCMC sampling
The mean, median, and mode are representative statistics of a distribution. In the case of a distribution of MCMC sampling, mean and median can be given by numpy functions mean() and percentile(). For the mode, we can estimate it by KDE estimation as follows
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]
An example of Python script, which estimate the mean, median, and mode (the mode here is called MAP (maximum a posteriori) estimate), is shown in Listing 1.
The function CalcMAPEst returns the mode, i.e. the MAP estimate for a sample from the posterior distribution.
The data is given as a list object in the script as follows.
data = [57, 30, 34, 59, 63, 26, 60,
58, 34, 56]
this data is analyzed by a univariate normal distribution model using the Stan script, which is given as a string object in the script as follows.
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);
}
}
"""
After MCMC sampling, the results are displayed as in Figure 1.
Figure 1
Means, medians, and modes are shown in the Figure.
import numpy as np
import pystan
import matplotlib.pyplot as plt
import seaborn as sb
import scipy.stats as ss
#
# The MAP
estimate, i.e.,the mode of a posterior distribution
#
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]
#
# Univariate
Normal Distribution Model
#
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)
fig, ax = plt.subplots(figsize = (12,
3))
mean_mu = fit['mu'].mean()
med_mu = np.percentile(fit['mu'], 50)
map_mu, max_v = CalcMAPEst(fit['mu'])
plt.subplot(1,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'Mode({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()
mean_sgm = fit['sgm'].mean()
med_sgm = np.percentile(fit['sgm'],
50)
map_sgm, max_v =
CalcMAPEst(fit['sgm'])
plt.subplot(1,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'Mode({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})')
plt.legend()
plt.tight_layout()
plt.savefig('Figure.png')
plt.show()