I would like to know if anyone has translated an example from Bayesian Data Analysis (3rd edition) to Stan. That example is Stratified sampling in pre-election polling, chapter 8, page 207. I know that there is a repository where you can find examples from BDA3 in Stan (here, but there are no examples from chapter 8) and I could do that, but I prefer to ask before I begin.
An example from BDA3 to Stan
Let’s begin with the nonhierarchical model. I will be using python. If you want to follow along with the code, first go here and download the file as .txt. I called it cbs_survey.txt.
Import some libraries:
import numpy as np
import matplotlib.pyplot as plt
import pystan
import pandas as pd
plt.style.use('seaborn-darkgrid')
plt.rc('font', size=12)
It is a survey of 1447 people, so:
participants = 1447
data = pd.read_csv('cbs_survey.txt', sep=' ', skiprows=2, skipinitialspace=True, index_col=False)
print(data)
We need the number of people of each region and each candidate.
data_obs = data[['bush', 'dukakis', 'other']].to_numpy()
print(data_obs)
And a variable called proportion:
proportion = data['proportion'].to_numpy() * participants
print(proportion)
And another variable:
valores = data_obs[:, :] * proportion.reshape(16, -1)
valores = np.round(valores)
np.sum(valores) # Check if the sum is equal to 1447
Now, Stan
model_non_hier="""
data {
int<lower=0> N;
int<lower=0> n;
int y_obs[N, n];
vector[n] alpha;
}
parameters {
simplex[n] theta[N];
}
model {
for (i in 1:N)
theta[i] ~ dirichlet(alpha);
for (i in 1:N)
y_obs[i] ~ multinomial(theta[i]);
}
"""
Now:
stan_model = pystan.StanModel(model_code=model_non_hier)
data2 = {'N': 16,
'n': 3,
'alpha': [1,1,1],
'y_obs': valores.astype(int)}
fit = stan_model.sampling(data=data2)
The output is very long.
print(fit.stansummary(digits_summary=5))
Inference for Stan model: anon_model_2a28c5090c429737936718a5aee45413.
4 chains, each with iter=2000; warmup=1000; thin=1;
post-warmup draws per chain=1000, total post-warmup draws=4000.
mean se_mean sd 2.5% 25% 50% 75% 97.5% n_eff Rhat
theta[1,1] 0.3007 0.0006 0.0643 0.1805 0.2566 0.2987 0.3413 0.4357 8933 0.9995
theta[2,1] 0.4902 0.0008 0.0726 0.3454 0.4426 0.4907 0.5382 0.6296 7457 0.9997
theta[3,1] 0.4648 0.0003 0.0379 0.3918 0.4387 0.4640 0.4902 0.5415 9417 0.9995
theta[4,1] 0.4587 0.0006 0.0566 0.3484 0.4198 0.4589 0.4974 0.5708 7798 0.9993
theta[5,1] 0.4011 0.0007 0.0690 0.2691 0.3541 0.3999 0.4479 0.5394 9671 0.9993
theta[6,1] 0.4424 0.0005 0.0510 0.3418 0.4083 0.4422 0.4758 0.5441 8823 0.9992
theta[7,1] 0.5041 0.0005 0.0467 0.4125 0.4723 0.5041 0.5361 0.5958 8730 0.9993
theta[8,1] 0.5474 0.0004 0.0412 0.4684 0.5191 0.5470 0.5772 0.6272 7904 0.9992
theta[9,1] 0.5420 0.0011 0.0984 0.3488 0.474 0.5443 0.6117 0.7277 7048 0.9994
theta[10,1] 0.4640 0.0005 0.0495 0.3675 0.4304 0.4641 0.4974 0.5617 8140 0.9994
theta[11,1] 0.5103 0.0006 0.0491 0.4110 0.4786 0.5098 0.5430 0.6038 6471 1.0000
theta[12,1] 0.5513 0.0004 0.0357 0.4815 0.5270 0.5506 0.5749 0.6215 7018 0.9993
theta[13,1] 0.4868 0.0008 0.0796 0.3304 0.4328 0.4865 0.5421 0.6381 8040 0.9992
theta[14,1] 0.5245 0.0006 0.0550 0.4146 0.4872 0.5255 0.5615 0.6307 7609 0.9998
theta[15,1] 0.5362 0.0004 0.0438 0.4496 0.5067 0.5373 0.5662 0.6197 9523 0.9996
theta[16,1] 0.5467 0.0005 0.0553 0.4364 0.5082 0.5471 0.5851 0.6519 8899 0.9997
theta[1,2] 0.5990 0.0008 0.0698 0.4569 0.5527 0.6008 0.6479 0.7298 7472 0.9993
theta[2,2] 0.4687 0.0008 0.0729 0.3284 0.4189 0.4679 0.5171 0.6174 7076 0.9997
theta[3,2] 0.4118 0.0004 0.0377 0.3386 0.3867 0.4115 0.4376 0.4858 8076 0.9997
theta[4,2] 0.5137 0.0006 0.0574 0.4014 0.4739 0.5140 0.5538 0.6282 7946 0.9994
theta[5,2] 0.4793 0.0007 0.0693 0.3462 0.4325 0.4794 0.5265 0.6144 8717 0.9993
theta[6,2] 0.4439 0.0005 0.0504 0.3453 0.4107 0.4431 0.4768 0.5447 8078 0.9993
theta[7,2] 0.3863 0.0004 0.0447 0.2973 0.3561 0.3854 0.4156 0.4760 8224 0.9993
theta[8,2] 0.3377 0.0004 0.0393 0.2636 0.3094 0.3373 0.3652 0.4153 7598 0.9993
theta[9,2] 0.2921 0.0010 0.0896 0.1357 0.2273 0.2865 0.3510 0.4831 7374 0.9993
theta[10,2] 0.4042 0.0005 0.0490 0.3086 0.3713 0.4038 0.4377 0.4997 7401 0.9997
theta[11,2] 0.4014 0.0006 0.0489 0.3059 0.3680 0.4007 0.4339 0.5005 6078 0.9996
theta[12,2] 0.3511 0.0004 0.0350 0.2808 0.3278 0.3510 0.3736 0.4225 7039 0.9996
theta[13,2] 0.4595 0.0009 0.0794 0.3089 0.4043 0.4587 0.5137 0.6184 7402 0.9993
theta[14,2] 0.3507 0.0006 0.0542 0.2478 0.3126 0.3489 0.3875 0.4623 8022 0.9994
theta[15,2] 0.3691 0.0004 0.0427 0.2854 0.3405 0.3685 0.3975 0.4549 7822 0.9995
theta[16,2] 0.3600 0.0005 0.0532 0.2594 0.3234 0.3589 0.3969 0.4663 8431 0.9996
theta[1,3] 0.1002 0.0004 0.0423 0.0333 0.0692 0.0952 0.1255 0.1975 7723 0.9993
theta[2,3] 0.0410 0.0003 0.0285 0.0051 0.0201 0.0347 0.0552 0.1128 6573 0.9999
theta[3,3] 0.1233 0.0002 0.0246 0.0793 0.1056 0.1223 0.1395 0.1750 8132 0.9999
theta[4,3] 0.0275 0.0002 0.0188 0.0032 0.0135 0.0236 0.0370 0.0749 6783 0.9997
theta[5,3] 0.1195 0.0005 0.0449 0.0459 0.0864 0.1152 0.1466 0.2226 8130 0.9995
theta[6,3] 0.1135 0.0003 0.0321 0.0577 0.0904 0.1112 0.1334 0.1824 7675 0.9992
theta[7,3] 0.1095 0.0003 0.0292 0.0603 0.0880 0.1073 0.1281 0.1701 7200 0.9994
theta[8,3] 0.1147 0.0003 0.0256 0.0689 0.0966 0.1132 0.1315 0.1691 7493 0.9994
theta[9,3] 0.1658 0.0009 0.0735 0.0503 0.1113 0.1582 0.2113 0.3322 6146 0.9995
theta[10,3] 0.1317 0.0004 0.0340 0.0724 0.1076 0.1297 0.1528 0.2042 7260 0.9997
theta[11,3] 0.0881 0.0003 0.0282 0.0409 0.0676 0.0855 0.1058 0.1517 7555 0.9993
theta[12,3] 0.0975 0.0002 0.0226 0.0583 0.0815 0.0959 0.1120 0.1456 6503 0.9997
theta[13,3] 0.0535 0.0004 0.0371 0.0060 0.0262 0.0458 0.0723 0.1464 6138 0.9994
theta[14,3] 0.1247 0.0004 0.0364 0.0619 0.0987 0.1217 0.148 0.2014 6789 0.9994
theta[15,3] 0.0946 0.0003 0.0266 0.0490 0.0754 0.0919 0.1113 0.1525 7251 0.9991
theta[16,3] 0.0932 0.0003 0.0315 0.0409 0.0706 0.0900 0.1123 0.1651 7889 0.9993
lp__ -1.4149e3 0.1127 4.2446-1.4241e3-1.4175e3-1.4145e3-1.4118e3-1.4076e3 1418 1.0038
Samples were drawn using NUTS at Tue Jul 30 17:33:35 2019.
For each parameter, n_eff is a crude measure of effective sample size,
and Rhat is the potential scale reduction factor on split chains (at
convergence, Rhat=1).
Let’s reproduce Figure 8.1(a).
samples = fit.extract(permuted=True)['theta']
print(sample.shape)
diff3 = []
for i in range(16):
result3 = samples[:, i, 0] - samples[:, i, 1]
diff3.append(list(result3))
diff3 = np.asarray(diff3)
res3 = np.sum(diff3.T * proportion / np.sum(proportion), axis=1)
plt.figure(figsize=(10, 6))
_, _, _ = plt.hist(res3, bins=20, edgecolor='w', density=True)
plt.show()
Nice!
Now the book goes to the hierarchical model (page 209).
As an aid to constructing a reasonable model, we transform the parameters to separate
partisan preference and probability of having a preference:\alpha_{1j} = \dfrac{\theta_{1j}}{\theta_{1j} + \theta_{2j}},
\alpha_{2j} = 1 - \theta_{3j}.
To get \alpha_{1j} and \alpha_{2j}:
alpha_2j = np.round((1 - data['other'].to_numpy()) * data.proportion.to_numpy() * participants)
print(alpha_2j)
alpha_1j = valores[:, 0] / (valores[:, 0] + valores[:, 1]) * data.proportion.to_numpy() * participants
print(alpha_1j)
new_values = np.round(np.stack([alpha_1j, alpha_2j], axis=1))
print(new_values)
With Stan:
model_hier="""
data {
int<lower=0> N;
int<lower=0> n;
int post[N, n];
}
parameters {
vector[n] mu;
cholesky_factor_corr[n] L;
simplex[n] alphas[N];
vector[n] beta[N];
}
transformed parameters{
vector[n] theta[N];
for (i in 1:N)
theta[i] = inv_logit(beta[i]);
}
model {
L ~ lkj_corr_cholesky(3.0);
mu ~ cauchy(0, 5);
beta ~ multi_normal_cholesky(mu, L);
for (i in 1:N)
alphas[i] ~ dirichlet(theta[i]);
for (i in 1:N)
post[i] ~ multinomial(alphas[i]);
}
generated quantities {
corr_matrix[n] Sigma;
Sigma = multiply_lower_tri_self_transpose(L);
}
"""
I have to say I am not happy about the code.
stan_model2 = pystan.StanModel(model_code=model_hier)
data3 = {'N': 16,
'n': 2,
'post': new_values.astype(int)}
fit2 = stan_model2.sampling(data=data3, algorithm='HMC', iter=4000, verbose=True)
Of course, there are warnings:
WARNING:pystan:n_eff / iter below 0.001 indicates that the effective sample size has likely been overestimated
WARNING:pystan:Rhat above 1.1 or below 0.9 indicates that the chains very likely have not mixed
WARNING:pystan:Skipping check of divergent transitions (divergence)
WARNING:pystan:Skipping check of transitions ending prematurely due to maximum tree depth limit (treedepth)
The output is very long
print(fit2.stansummary(digits_summary=5))
Inference for Stan model: anon_model_958bb8552d6652a8814e6044bd42f4a8.
4 chains, each with iter=4000; warmup=2000; thin=1;
post-warmup draws per chain=2000, total post-warmup draws=8000.
mean se_mean sd 2.5% 25% 50% 75% 97.5% n_eff Rhat
mu[1] 2.2705e1 8.01884.6381e1 1.1962 3.6407 7.20241.7178e11.7396e2 33 1.1352
mu[2] 2.0993e1 4.87353.433e1 2.2160 4.9229 8.72791.931e11.4692e2 50 1.0651
L[1,1] 1.0 nan 0.0 1.0 1.0 1.0 1.0 1.0 nan nan
L[2,1] -0.000 0.0056 0.3769 -0.709 -0.275 0.0006 0.2848 0.7012 4423 0.9999
L[1,2] 0.0 nan 0.0 0.0 0.0 0.0 0.0 0.0 nan nan
L[2,2] 0.9210 0.0023 0.0976 0.6453 0.8878 0.9600 0.9904 0.9999 1741 1.0020
alphas[1,1] 0.2707 0.0002 0.0572 0.1649 0.2311 0.2686 0.3065 0.3915 39793 0.9996
alphas[2,1] 0.352 0.0006 0.0556 0.2468 0.3122 0.3519 0.3904 0.4620 8197 0.9998
alphas[3,1] 0.3770 0.0001 0.0308 0.3166 0.3565 0.3769 0.3975 0.4385 34281 0.9995
alphas[4,1] 0.3306 0.0005 0.0456 0.2474 0.2974 0.3295 0.3627 0.4174 7888 1.0011
alphas[5,1] 0.3426 0.0010 0.0577 0.2351 0.3024 0.3405 0.3821 0.4602 2831 1.0003
alphas[6,1] 0.3604 0.0002 0.0361 0.2878 0.3384 0.3600 0.3816 0.4350 19799 0.9995
alphas[7,1] 0.3890 0.0003 0.0362 0.3204 0.3638 0.3891 0.4134 0.4601 13714 0.9997
alphas[8,1] 0.4110 0.0005 0.0340 0.3451 0.3896 0.4106 0.4324 0.48 4225 1.0001
alphas[9,1] 0.4282 0.0004 0.0822 0.2730 0.3697 0.4263 0.4849 0.5934 33091 0.9995
alphas[10,1] 0.3791 0.0009 0.0393 0.3048 0.3516 0.3782 0.4053 0.4580 1595 1.0016
alphas[11,1] 0.3804 0.0003 0.0379 0.3069 0.3543 0.3802 0.4063 0.4566 9941 0.9996
alphas[12,1] 0.4053 0.0002 0.0291 0.3484 0.3852 0.4051 0.425 0.4637 22000 0.9997
alphas[13,1] 0.3529 0.0004 0.0701 0.2175 0.3060 0.3507 0.3983 0.4993 28748 0.9996
alphas[14,1] 0.4052 0.0002 0.0456 0.3173 0.3727 0.4055 0.4370 0.4955 33739 0.9995
alphas[15,1] 0.3970 0.0001 0.0325 0.3342 0.3747 0.3964 0.4189 0.4611 49240 0.9996
alphas[16,1] 0.3979 0.0004 0.0453 0.3108 0.3667 0.3974 0.4286 0.4883 11043 0.9997
alphas[1,2] 0.7292 0.0002 0.0572 0.6084 0.6934 0.7313 0.7688 0.8350 39793 0.9996
alphas[2,2] 0.648 0.0006 0.0556 0.5379 0.6095 0.6480 0.6877 0.7531 8197 0.9998
alphas[3,2] 0.6229 0.0001 0.0308 0.5614 0.6024 0.6230 0.6434 0.6833 34281 0.9995
alphas[4,2] 0.6693 0.0005 0.0456 0.5825 0.6372 0.6704 0.7025 0.7525 7888 1.0011
alphas[5,2] 0.6573 0.0010 0.0577 0.5397 0.6179 0.6594 0.6975 0.7648 2831 1.0003
alphas[6,2] 0.6395 0.0002 0.0361 0.5649 0.6183 0.6399 0.6615 0.7121 19799 0.9995
alphas[7,2] 0.6109 0.0003 0.0362 0.5398 0.5865 0.6108 0.6362 0.6796 13714 0.9997
alphas[8,2] 0.5889 0.0005 0.0340 0.52 0.5675 0.5893 0.6104 0.6548 4225 1.0001
alphas[9,2] 0.5717 0.0004 0.0822 0.4065 0.5150 0.5736 0.6302 0.7269 33091 0.9995
alphas[10,2] 0.6208 0.0009 0.0393 0.5419 0.5946 0.6217 0.6483 0.6951 1595 1.0016
alphas[11,2] 0.6195 0.0003 0.0379 0.5433 0.5937 0.6197 0.6456 0.6930 9941 0.9996
alphas[12,2] 0.5947 0.0002 0.0291 0.5362 0.575 0.5948 0.6147 0.6515 22000 0.9997
alphas[13,2] 0.6470 0.0004 0.0701 0.5006 0.6016 0.6493 0.6939 0.7824 28748 0.9996
alphas[14,2] 0.5948 0.0002 0.0456 0.5044 0.5629 0.5944 0.6272 0.6827 33739 0.9995
alphas[15,2] 0.6029 0.0001 0.0325 0.5388 0.5810 0.6035 0.6253 0.6657 49240 0.9996
alphas[16,2] 0.6020 0.0004 0.0453 0.5116 0.5713 0.6025 0.6333 0.6891 11043 0.9997
beta[1,1] 2.2713e1 8.01864.6367e1 0.7253 3.7020 7.29131.7041e11.7394e2 33 1.1354
beta[2,1] 2.2719e1 8.01924.6396e1 0.7869 3.6451 7.28721.7319e11.7264e2 33 1.1352
beta[3,1] 2.2694e1 8.02054.6386e1 0.7152 3.6263 7.24291.7117e11.7349e2 33 1.1352
beta[4,1] 2.2711e1 8.01274.6375e1 0.7990 3.7203 7.2744 1.72e11.7371e2 33 1.1350
beta[5,1] 2.2706e1 8.01994.6389e1 0.6911 3.6927 7.28541.7143e11.7308e2 33 1.1352
beta[6,1] 2.2696e1 8.01374.6386e1 0.7916 3.6592 7.29001.7171e11.7384e2 34 1.1350
beta[7,1] 2.2715e1 8.02274.6393e1 0.7666 3.7174 7.29451.7197e11.7365e2 33 1.1354
beta[8,1] 2.2714e1 8.01464.6374e1 0.7420 3.7158 7.28731.7187e11.7354e2 33 1.1352
beta[9,1] 2.2729e1 8.01774.6387e1 0.8191 3.6922 7.16941.7145e11.7301e2 33 1.1352
beta[10,1] 2.2718e1 8.01684.6386e1 0.8272 3.7168 7.22621.7146e11.744e2 33 1.1351
beta[11,1] 2.2725e1 8.01984.6402e1 0.8458 3.6746 7.2161.7197e11.7302e2 33 1.1351
beta[12,1] 2.2711e1 8.01864.638e1 0.6594 3.6695 7.27421.7207e11.7349e2 33 1.1352
beta[13,1] 2.2708e1 8.01604.637e1 0.6844 3.7449 7.25291.7199e11.7331e2 33 1.1352
beta[14,1] 2.2719e1 8.01824.6392e1 0.7666 3.6791 7.27351.7018e11.735e2 33 1.1352
beta[15,1] 2.2716e1 8.01634.6394e1 0.8257 3.6758 7.26291.7119e11.7331e2 33 1.1351
beta[16,1] 2.2715e1 8.01634.6387e1 0.7689 3.6845 7.22611.7126e11.731e2 33 1.1350
beta[1,2] 2.101e1 4.87323.4339e1 1.8136 4.9731 8.79371.943e11.4617e2 50 1.0649
beta[2,2] 2.1011e1 4.87093.4341e1 1.7563 5.0524 8.80251.9373e11.4657e2 50 1.0651
beta[3,2] 2.0995e1 4.87163.4329e1 1.6634 5.0308 8.7981.9455e11.4667e2 50 1.0649
beta[4,2] 2.1008e1 4.87323.4328e1 1.8091 5.0136 8.75131.9424e11.4693e2 50 1.0652
beta[5,2] 2.1034e1 4.86883.4318e1 1.8465 5.0907 8.73281.9448e11.4711e2 50 1.0650
beta[6,2] 2.1014e1 4.87743.4345e1 1.7678 5.0475 8.74091.9357e11.4671e2 50 1.0652
beta[7,2] 2.1019e1 4.87293.4342e1 1.8202 5.0471 8.77261.9426e11.4751e2 50 1.0650
beta[8,2] 2.1e1 4.87203.4343e1 1.7964 4.9918 8.78851.937e11.4724e2 50 1.0648
beta[9,2] 2.1003e1 4.87543.434e1 1.7889 4.9947 8.77091.9506e11.467e2 50 1.0651
beta[10,2] 2.1019e1 4.87453.4355e1 1.7801 5.0105 8.83531.9436e11.4755e2 50 1.0651
beta[11,2] 2.1011e1 4.88323.436e1 1.7745 4.9971 8.81831.9411e11.4697e2 50 1.0652
beta[12,2] 2.1e1 4.86893.4323e1 1.7969 5.0186 8.76581.9163e11.4668e2 50 1.0650
beta[13,2] 2.1024e1 4.86963.4322e1 1.7591 5.0657 8.77951.9546e11.4676e2 50 1.0645
beta[14,2] 2.1002e1 4.87323.4324e1 1.8039 4.9603 8.76831.9452e11.4686e2 50 1.0652
beta[15,2] 2.1013e1 4.87403.434e1 1.8079 4.9877 8.83411.9378e11.466e2 50 1.0654
beta[16,2] 2.1012e1 4.87293.4341e1 1.7040 5.0002 8.79891.9428e11.4713e2 50 1.0648
theta[1,1] 0.9628 0.0015 0.0904 0.6737 0.9759 0.9993 1.0 1.0 3549 1.0039
theta[2,1] 0.9627 0.0016 0.0897 0.6871 0.9745 0.9993 1.0 1.0 3053 1.0045
theta[3,1] 0.9631 0.0016 0.0900 0.6715 0.9740 0.9992 1.0 1.0 2994 1.0047
theta[4,1] 0.9642 0.0013 0.0869 0.6897 0.9763 0.9993 1.0 1.0 4015 1.0027
theta[5,1] 0.9629 0.0016 0.0884 0.6662 0.9757 0.9993 1.0 1.0 2989 1.0044
theta[6,1] 0.9628 0.0013 0.0896 0.6881 0.9749 0.9993 1.0 1.0 4184 1.0038
theta[7,1] 0.9633 0.0014 0.0885 0.6827 0.9762 0.9993 1.0 1.0 3702 1.0046
theta[8,1] 0.9629 0.0013 0.0904 0.6774 0.9762 0.9993 1.0 1.0 4262 1.0028
theta[9,1] 0.9652 0.0013 0.0832 0.6940 0.9756 0.9992 1.0 1.0 3682 1.0043
theta[10,1] 0.9645 0.0014 0.0858 0.6957 0.9762 0.9992 1.0 1.0 3582 1.0043
theta[11,1] 0.9650 0.0012 0.0827 0.6997 0.9752 0.9992 1.0 1.0 4558 1.0042
theta[12,1] 0.9627 0.0015 0.0897 0.6591 0.9751 0.9993 1.0 1.0 3395 1.0047
theta[13,1] 0.9627 0.0014 0.0899 0.6647 0.9769 0.9992 1.0 1.0 3858 1.0052
theta[14,1] 0.9632 0.0013 0.0893 0.6827 0.9753 0.9993 1.0 1.0 4132 1.0028
theta[15,1] 0.9643 0.0013 0.0858 0.6954 0.9753 0.9993 1.0 1.0 3811 1.0041
theta[16,1] 0.9638 0.0015 0.0868 0.6833 0.9755 0.9992 1.0 1.0 3231 1.0035
theta[1,2] 0.9855 0.0006 0.0429 0.8598 0.9931 0.9998 1.0 1.0 4995 1.0018
theta[2,2] 0.9848 0.0007 0.0458 0.8527 0.9936 0.9998 1.0 1.0 4314 1.0017
theta[3,2] 0.9846 0.0007 0.0454 0.8407 0.9935 0.9998 1.0 1.0 3969 1.0011
theta[4,2] 0.985 0.0006 0.0472 0.8592 0.9934 0.9998 1.0 1.0 5750 1.0013
theta[5,2] 0.9852 0.0007 0.0464 0.8637 0.9938 0.9998 1.0 1.0 4179 1.0008
theta[6,2] 0.9850 0.0006 0.0466 0.8541 0.9936 0.9998 1.0 1.0 5104 1.0019
theta[7,2] 0.9856 0.0006 0.0434 0.8606 0.9936 0.9998 1.0 1.0 4667 1.0007
theta[8,2] 0.9851 0.0006 0.0454 0.8577 0.9932 0.9998 1.0 1.0 4716 1.0008
theta[9,2] 0.9848 0.0006 0.0482 0.8567 0.9932 0.9998 1.0 1.0 5013 1.0013
theta[10,2] 0.9848 0.0006 0.0474 0.8557 0.9933 0.9998 1.0 1.0 5051 1.0013
theta[11,2] 0.9853 0.0006 0.0439 0.8550 0.9932 0.9998 1.0 1.0 4887 1.0009
theta[12,2] 0.9851 0.0006 0.0463 0.8577 0.9934 0.9998 1.0 1.0 5339 1.0006
theta[13,2] 0.9851 0.0006 0.0453 0.8531 0.9937 0.9998 1.0 1.0 5175 1.0003
theta[14,2] 0.9855 0.0006 0.0441 0.8586 0.9930 0.9998 1.0 1.0 4960 1.0016
theta[15,2] 0.9851 0.0006 0.0450 0.8591 0.9932 0.9998 1.0 1.0 4773 1.0011
theta[16,2] 0.9848 0.0007 0.0466 0.8460 0.9933 0.9998 1.0 1.0 4484 1.0013
Sigma[1,1] 1.0 nan 0.0 1.0 1.0 1.0 1.0 1.0 nan nan
Sigma[2,1] -0.000 0.0056 0.3769 -0.709 -0.275 0.0006 0.2848 0.7012 4423 0.9999
Sigma[1,2] -0.000 0.0056 0.3769 -0.709 -0.275 0.0006 0.2848 0.7012 4423 0.9999
Sigma[2,2] 1.01.0388e-188.708e-17 1.0 1.0 1.0 1.0 1.0 7027 0.9995
lp__ -1.4472e3 0.2517 5.6630-1.4594e3-1.4508e3-1.4469e3-1.4434e3-1.437e3 506 1.0095
Samples were drawn using HMC at Tue Jul 30 20:01:26 2019.
For each parameter, n_eff is a crude measure of effective sample size,
and Rhat is the potential scale reduction factor on split chains (at
convergence, Rhat=1).
Let’s reproduce Figure 8.1(b)
samples2 = fit2.extract(permuted=True)['alphas']
th5 = []
for i in range(16):
result5 = 2 * samples2[:, i, 0] * samples2[:, i, 1] - samples2[:, i, 1]
th5.append(list(result5))
th5 = np.asarray(th5)
print(th5.shape)
res5 = np.sum(th5.T * proportion / np.sum(proportion), axis=1)
plt.figure(figsize=(10, 6))
_, _, _ = plt.hist(res5 , bins=20, edgecolor='w', density=True)
plt.show()
Do you know why that is weird?
The correct figure is (I think):
I have only changed the sign of res5. My Stan code could be wrong… but when you get the same figure with pymc3, you begin to doubt. The goal is then to improve the model (the diagnostics are not good), and that will be very hard.