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.
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.