Comparison of Point Estimates in Bayesian Analyisi
Means, medians, and modes of posterior distributions
Yasuharu Okamoto
When we want to estimate values for parameters in a model, there are several methods. Methods of moments or maximum likelihood estimates (MLEs) are popular ones. If we use a Bayesian method, we must estimate the value for a parameter from the posterior distribution. Means, medians, and modes are popular statistics, which represent the position of the distribution. Different researches recommend different statistics. Means are recommended by Lambert (2018). But, Gelman et al. (2021) criticized means as unstable, and recommend medians. Furthermore, Gelman et al. (2021) say “the posterior mode, which provides the best overall fit to data and prior distribution (p. 107).”
Okamoto (2022) compared the three statistics and concluded that the mode is best, as shown in the following.
Okamoto (in press) discussed the three models, a binomial model, a univariate normal distribution, and a bivariate normal distribution. The script files, which were edited versions of the scripts used in the calculations for Okamoto (in press), were archived in the file PointEstimates.zip. This file PointEstimates.zip can be freely downloaded and used on user’s responsibility
.
In the case of a binomial model,
the posterior distribution for the uniform prior distribution is obtained as a beta distribution. As shown in Figure 1, the mode seems to be the best.
Figure 1
The graph in Figure 1 was drawn by the Python script in Listing 1.
To compare the three estimates, biases and
RMSEs (root mean square error) were calculated for sample sizes N=100 and
N=200. Figures 2, 3 and 4 show the results for true value 
. These figures were drawn by the script in Listing
2.
Figure 2
Figure 3
Figure 4
In the case of a normal distribution, the following model was set:
Run the script in Listing 3, we get the results in the file ResultsCheck.txt as follows:
*/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
Actual results vary from the results shown in Table 1 of Okamoto(in press). The differences between these values and the values in Table 1 of Okamoto (in press) can be accepted when considering the values SE.
In the case of a bivariate normal distribution, the following model was set:
When we run the script in Listing 5, we get the results in the file ResultsN50.txt like the following:
*/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
Values in each line denote the values for
parameter values 
to 0.9. These values are
not exactly the same as in Tables 2 and 3. These discrepancies are due to the
method of the simulation, in which random number generator is used, so the
results vary from simulation to simulation. The values of SEs show that the
discrepancies between the above values and values in Tables 2 and 3 are within
variations due to a random number generator.
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. (2022). Some comparisons of point estimates of parameters by Bayesian analysis with respect to RMSE and bias: A report from the point of view of experimental psychology. The Japanese Journal of Psychonomic Science, vol. 41, No. 1. https://doi.org/10.14947/psychono.41.1
Listing 1. The script for drawing the graph in Figure 1.
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()
Listing 2. The script for Figures 2, 3 and 4.
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()
Listing 3. The script for simulation of a univariate normal distribution, which uses Stan script in Listing 4.
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()
Listing 4. A Stan script for a univariate normal distribution.
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);
}
}
Listing 5. Simulation for a bivariate
normal distribution, which uses the Stan script in Listing 6.
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()
Listing 6. A Stan script for a bivariate normal distribution
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);
}
}