import arviz as az
import matplotlib.pyplot as plt
import numpy as np
import pandas as pd
import pymc3 as pm
import seaborn as sns
import theano.tensor as tt
import warnings
from scipy import stats
from scipy.special import expit as logistic
%config InlineBackend.figure_format = 'retina'
az.style.use('arviz-darkgrid')
warnings.simplefilter(action="ignore", category=FutureWarning)
d = pd.read_csv('Data/chimpanzees.csv', sep=";")
# we change "actor" to zero-index
d.actor = d.actor - 1
with pm.Model() as model_10_1:
a = pm.Normal('a', 0, 10)
bp = pm.Normal('bp', 0, 10)
p = pm.math.invlogit(a)
pulled_left = pm.Binomial('pulled_left', 1, p, observed=d.pulled_left)
trace_10_1 = pm.sample(1000, tune=1000)
Auto-assigning NUTS sampler... Initializing NUTS using jitter+adapt_diag... Multiprocess sampling (2 chains in 2 jobs) NUTS: [bp, a]
HBox(children=(IntProgress(value=0, description='Sampling 2 chains', max=4000, style=ProgressStyle(description…
df_10_1 = az.summary(trace_10_1, credible_interval=.89, round_to=2)
df_10_1
mean | sd | hpd_5.5% | hpd_94.5% | mcse_mean | mcse_sd | ess_mean | ess_sd | ess_bulk | ess_tail | r_hat | |
---|---|---|---|---|---|---|---|---|---|---|---|
a | 0.32 | 0.09 | 0.18 | 0.47 | 0.00 | 0.00 | 2195.32 | 2047.05 | 2215.79 | 1408.13 | 1.0 |
bp | 0.02 | 9.97 | -14.54 | 16.74 | 0.21 | 0.21 | 2186.02 | 1134.28 | 2190.82 | 1583.60 | 1.0 |
logistic(df_10_1.iloc[:,2:4]).round(5)
hpd_5.5% | hpd_94.5% | |
---|---|---|
a | 0.54488 | 0.61538 |
bp | 0.00000 | 1.00000 |
with pm.Model() as model_10_2:
a = pm.Normal('a', 0, 10)
bp = pm.Normal('bp', 0, 10)
p = pm.math.invlogit(a + bp * d.prosoc_left)
pulled_left = pm.Binomial('pulled_left', 1, p, observed=d.pulled_left)
trace_10_2 = pm.sample(1000, tune=1000)
with pm.Model() as model_10_3:
a = pm.Normal('a', 0, 10)
bp = pm.Normal('bp', 0, 10)
bpC = pm.Normal('bpC', 0, 10)
p = pm.math.invlogit(a + (bp + bpC * d.condition) * d.prosoc_left)
pulled_left = pm.Binomial('pulled_left', 1, p, observed=d.pulled_left)
trace_10_3 = pm.sample(1000, tune=1000)
Auto-assigning NUTS sampler... Initializing NUTS using jitter+adapt_diag... Multiprocess sampling (2 chains in 2 jobs) NUTS: [bp, a]
HBox(children=(IntProgress(value=0, description='Sampling 2 chains', max=4000, style=ProgressStyle(description…
Auto-assigning NUTS sampler... Initializing NUTS using jitter+adapt_diag... Multiprocess sampling (2 chains in 2 jobs) NUTS: [bpC, bp, a]
HBox(children=(IntProgress(value=0, description='Sampling 2 chains', max=4000, style=ProgressStyle(description…
comp_df = az.compare({'m10.1' : trace_10_1,
'm10.2' : trace_10_2,
'm10.3' : trace_10_3})
comp_df
waic | p_waic | d_waic | weight | se | dse | warning | waic_scale | |
---|---|---|---|---|---|---|---|---|
m10.2 | 680.338 | 1.91809 | 0 | 0.635535 | 9.17014 | 0 | False | deviance |
m10.3 | 682.377 | 3.01781 | 2.0389 | 0.239546 | 7.01294 | 0.848115 | False | deviance |
m10.1 | 687.922 | 0.990988 | 7.58409 | 0.124919 | 9.15942 | 6.32352 | False | deviance |
az.plot_compare(comp_df);
az.summary(trace_10_3, credible_interval=.89, round_to=2)
mean | sd | hpd_5.5% | hpd_94.5% | mcse_mean | mcse_sd | ess_mean | ess_sd | ess_bulk | ess_tail | r_hat | |
---|---|---|---|---|---|---|---|---|---|---|---|
a | 0.04 | 0.13 | -0.15 | 0.25 | 0.00 | 0.00 | 1224.98 | 1073.93 | 1222.69 | 1338.51 | 1.0 |
bp | 0.62 | 0.22 | 0.28 | 0.98 | 0.01 | 0.00 | 1081.41 | 1040.40 | 1084.76 | 1328.61 | 1.0 |
bpC | -0.11 | 0.26 | -0.54 | 0.27 | 0.01 | 0.01 | 1108.54 | 999.90 | 1103.29 | 1371.47 | 1.0 |
np.exp(0.61)
1.8404313987816374
logistic(4)
0.9820137900379085
logistic(4 + 0.61)
0.9901462444767687
d_pred = pd.DataFrame({'prosoc_left' : [0, 1, 0, 1], 'condition' : [0, 0, 1, 1]})
traces = [trace_10_1, trace_10_2, trace_10_3]
models = [model_10_1, model_10_2, model_10_3]
chimp_ensemble = pm.sample_ppc_w(traces=traces, models=models, samples=1000, weights=comp_df.weight.sort_index(ascending=True))
/home/osvaldo/anaconda3/lib/python3.7/site-packages/ipykernel_launcher.py:6: DeprecationWarning: sample_ppc() is deprecated. Please use sample_posterior_predictive_w()
HBox(children=(IntProgress(value=0, max=1000), HTML(value='')))
rt = chimp_ensemble['pulled_left']
pred_mean = np.zeros((1000, 4))
cond = d.condition.unique()
prosoc_l = d.prosoc_left.unique()
for i in range(len(rt)):
tmp = []
if rt[i].size < 2:
continue
for cp in cond:
for pl in prosoc_l:
tmp.append(np.mean(rt[i][(d.prosoc_left==pl) & (d.chose_prosoc==cp)]))
pred_mean[i] = tmp
ticks = range(4)
mp = pred_mean.mean(0)
az.plot_hpd(ticks, pred_mean, color='k', smooth=False)
plt.plot(mp, color='k')
plt.xticks(ticks, ("0/0","1/0","0/1","1/1"))
chimps = d.groupby(['actor', 'prosoc_left', 'condition']).agg('mean')['pulled_left'].values.reshape(7, -1)
for i in range(7):
plt.plot(chimps[i], 'C0')
plt.ylim(0, 1.1);
This is the same as 10.6, but in the book using MCMC rather than quadratic approximation.
with pm.Model() as model_10_4:
a = pm.Normal('alpha', 0, 10, shape=len(d.actor.unique()))
bp = pm.Normal('bp', 0, 10)
bpC = pm.Normal('bpC', 0, 10)
p = pm.math.invlogit(a[d.actor.values] + (bp + bpC * d.condition) * d.prosoc_left)
pulled_left = pm.Binomial('pulled_left', 1, p, observed=d.pulled_left)
trace_10_4 = pm.sample(1000, tune=1000)
Auto-assigning NUTS sampler... Initializing NUTS using jitter+adapt_diag... Multiprocess sampling (2 chains in 2 jobs) NUTS: [bpC, bp, alpha]
HBox(children=(IntProgress(value=0, description='Sampling 2 chains', max=4000, style=ProgressStyle(description…
There were 3 divergences after tuning. Increase `target_accept` or reparameterize. There were 10 divergences after tuning. Increase `target_accept` or reparameterize.
# remember we use a zero-index
d['actor'].unique()
array([0, 1, 2, 3, 4, 5, 6])
az.summary(trace_10_4, credible_interval=.89, round_to=2)
mean | sd | hpd_5.5% | hpd_94.5% | mcse_mean | mcse_sd | ess_mean | ess_sd | ess_bulk | ess_tail | r_hat | |
---|---|---|---|---|---|---|---|---|---|---|---|
alpha[0] | -0.73 | 0.29 | -1.19 | -0.29 | 0.01 | 0.01 | 1504.52 | 1412.24 | 1507.14 | 1342.57 | 1.0 |
alpha[1] | 10.65 | 5.26 | 3.37 | 17.68 | 0.19 | 0.15 | 741.37 | 600.19 | 1093.08 | 790.00 | 1.0 |
alpha[2] | -1.04 | 0.29 | -1.48 | -0.58 | 0.01 | 0.01 | 1233.41 | 1233.41 | 1225.45 | 1308.71 | 1.0 |
alpha[3] | -1.05 | 0.29 | -1.54 | -0.61 | 0.01 | 0.01 | 1359.02 | 1289.46 | 1366.34 | 1260.81 | 1.0 |
alpha[4] | -0.74 | 0.27 | -1.17 | -0.30 | 0.01 | 0.01 | 1501.36 | 1416.50 | 1505.25 | 1523.97 | 1.0 |
alpha[5] | 0.23 | 0.27 | -0.16 | 0.69 | 0.01 | 0.01 | 1622.82 | 1165.74 | 1631.27 | 1218.23 | 1.0 |
alpha[6] | 1.80 | 0.39 | 1.23 | 2.46 | 0.01 | 0.01 | 1603.79 | 1397.90 | 1691.53 | 1193.19 | 1.0 |
bp | 0.84 | 0.27 | 0.41 | 1.27 | 0.01 | 0.01 | 1009.23 | 967.62 | 1012.28 | 1118.15 | 1.0 |
bpC | -0.15 | 0.30 | -0.64 | 0.31 | 0.01 | 0.01 | 1470.11 | 1148.59 | 1470.57 | 1309.26 | 1.0 |
post = pm.trace_to_dataframe(trace_10_4)
post.head()
alpha__0 | alpha__1 | alpha__2 | alpha__3 | alpha__4 | alpha__5 | alpha__6 | bp | bpC | |
---|---|---|---|---|---|---|---|---|---|
0 | -0.968013 | 18.225467 | -1.200813 | -1.161509 | -0.410280 | 0.188969 | 1.105035 | 0.881806 | -0.279611 |
1 | -0.931144 | 9.180765 | -1.187670 | -1.205491 | -1.529347 | 0.006999 | 2.097104 | 1.172396 | -0.389430 |
2 | -0.792369 | 5.581284 | -0.489599 | -1.386701 | -1.135585 | 0.220857 | 1.689971 | 0.949201 | 0.038198 |
3 | -0.860857 | 6.848895 | -1.672781 | -1.017318 | -0.549615 | 0.397582 | 1.740100 | 0.814674 | -0.444185 |
4 | -0.904558 | 14.308174 | -1.148785 | -1.330190 | -1.321857 | -0.126976 | 1.631864 | 1.390591 | -0.445198 |
az.plot_kde(post['alpha__1']);
rt = pm.sample_posterior_predictive(trace_10_4, 1000, model_10_4)['pulled_left']
/home/osvaldo/proyectos/00_PyMC3/pymc3/pymc3/sampling.py:1100: UserWarning: samples parameter is smaller than nchains times ndraws, some draws and/or chains may not be represented in the returned posterior predictive sample warnings.warn("samples parameter is smaller than nchains times ndraws, some draws "
HBox(children=(IntProgress(value=0, max=1000), HTML(value='')))
chimp = 2
pred_mean = np.zeros((1000, 4))
cond = d.condition.unique()
prosoc_l = d.prosoc_left.unique()
for i in range(len(rt)):
tmp = []
for cp in cond:
for pl in prosoc_l:
tmp.append(np.mean(rt[i][(d.prosoc_left == pl) &
(d.chose_prosoc == cp) &
(d.actor == chimp)]))
pred_mean[i] = tmp
ticks = range(4)
mp = pred_mean.mean(0)
hpd = pm.hpd(pred_mean, alpha=0.11)
plt.fill_between(ticks, hpd[:,1], hpd[:,0], alpha=0.25, color='k')
plt.plot(mp, color='k')
plt.xticks(ticks, ("0/0","1/0","0/1","1/1"))
chimps = d[d.actor == chimp].groupby(['condition', 'prosoc_left', ]).agg('mean')['pulled_left'].values
plt.plot(chimps, 'C0')
plt.ylim(0, 1.1);
d_aggregated = d.groupby(['actor', 'condition', 'prosoc_left', ])['pulled_left'].sum().reset_index()
d_aggregated.head(8)
actor | condition | prosoc_left | pulled_left | |
---|---|---|---|---|
0 | 0 | 0 | 0 | 6 |
1 | 0 | 0 | 1 | 9 |
2 | 0 | 1 | 0 | 5 |
3 | 0 | 1 | 1 | 10 |
4 | 1 | 0 | 0 | 18 |
5 | 1 | 0 | 1 | 18 |
6 | 1 | 1 | 0 | 18 |
7 | 1 | 1 | 1 | 18 |
with pm.Model() as model_10_5:
a = pm.Normal('alpha', 0, 10)
bp = pm.Normal('bp', 0, 10)
bpC = pm.Normal('bpC', 0, 10)
p = pm.math.invlogit(a + (bp + bpC * d_aggregated.condition) * d_aggregated.prosoc_left)
pulled_left = pm.Binomial('pulled_left', 18, p, observed=d_aggregated.pulled_left)
trace_10_5 = pm.sample(1000, tune=1000)
Auto-assigning NUTS sampler... Initializing NUTS using jitter+adapt_diag... Multiprocess sampling (2 chains in 2 jobs) NUTS: [bpC, bp, alpha]
HBox(children=(IntProgress(value=0, description='Sampling 2 chains', max=4000, style=ProgressStyle(description…
az.summary(trace_10_5, credible_interval=.89, round_to=2)
mean | sd | hpd_5.5% | hpd_94.5% | mcse_mean | mcse_sd | ess_mean | ess_sd | ess_bulk | ess_tail | r_hat | |
---|---|---|---|---|---|---|---|---|---|---|---|
alpha | 0.05 | 0.12 | -0.15 | 0.24 | 0.00 | 0.00 | 1382.46 | 907.65 | 1383.81 | 1072.35 | 1.0 |
bp | 0.61 | 0.23 | 0.22 | 0.96 | 0.01 | 0.01 | 973.68 | 973.68 | 979.60 | 803.54 | 1.0 |
bpC | -0.10 | 0.27 | -0.53 | 0.35 | 0.01 | 0.01 | 1087.00 | 877.87 | 1085.07 | 1008.47 | 1.0 |
# hacky check of similarity to 10_3, within a hundreth
np.isclose(az.summary(trace_10_5), az.summary(trace_10_3), atol=0.01)
array([[ True, True, False, False, True, True, False, False, False, False, True], [False, True, False, False, True, True, False, False, False, False, True], [False, False, False, False, True, True, False, False, False, False, True]])
d_ad = pd.read_csv('./Data/UCBadmit.csv', sep=';')
d_ad.head(8)
dept | applicant.gender | admit | reject | applications | |
---|---|---|---|---|---|
1 | A | male | 512 | 313 | 825 |
2 | A | female | 89 | 19 | 108 |
3 | B | male | 353 | 207 | 560 |
4 | B | female | 17 | 8 | 25 |
5 | C | male | 120 | 205 | 325 |
6 | C | female | 202 | 391 | 593 |
7 | D | male | 138 | 279 | 417 |
8 | D | female | 131 | 244 | 375 |
d_ad['male'] = (d_ad['applicant.gender'] == 'male').astype(int)
with pm.Model() as model_10_6:
a = pm.Normal('a', 0, 10)
bm = pm.Normal('bm', 0, 10)
p = pm.math.invlogit(a + bm * d_ad.male)
admit = pm.Binomial('admit', p=p, n=d_ad.applications, observed=d_ad.admit)
trace_10_6 = pm.sample(1000, tune=1000)
with pm.Model() as model_10_7:
a = pm.Normal('a', 0, 10)
p = pm.math.invlogit(a)
admit = pm.Binomial('admit', p=p, n=d_ad.applications, observed=d_ad.admit)
trace_10_7 = pm.sample(1000, tune=1000)
Auto-assigning NUTS sampler... Initializing NUTS using jitter+adapt_diag... Multiprocess sampling (2 chains in 2 jobs) NUTS: [bm, a]
HBox(children=(IntProgress(value=0, description='Sampling 2 chains', max=4000, style=ProgressStyle(description…
Auto-assigning NUTS sampler... Initializing NUTS using jitter+adapt_diag... Multiprocess sampling (2 chains in 2 jobs) NUTS: [a]
HBox(children=(IntProgress(value=0, description='Sampling 2 chains', max=4000, style=ProgressStyle(description…
# Something goofy here...
# not even close to WAIC values, larger standard error
comp_df = az.compare({'m10.6' : trace_10_6,
'm10.7' : trace_10_7})
comp_df
/home/osvaldo/proyectos/00_PyMC3/arviz/arviz/stats/stats.py:1078: UserWarning: For one or more samples the posterior variance of the log predictive densities exceeds 0.4. This could be indication of WAIC starting to fail see http://arxiv.org/abs/1507.04544 for details """ /home/osvaldo/proyectos/00_PyMC3/arviz/arviz/stats/stats.py:1078: UserWarning: For one or more samples the posterior variance of the log predictive densities exceeds 0.4. This could be indication of WAIC starting to fail see http://arxiv.org/abs/1507.04544 for details """
waic | p_waic | d_waic | weight | se | dse | warning | waic_scale | |
---|---|---|---|---|---|---|---|---|
m10.6 | 990.347 | 111.824 | 0 | 0.614148 | 316.551 | 0 | True | deviance |
m10.7 | 1041.33 | 79.544 | 50.9798 | 0.385852 | 301.674 | 157.884 | True | deviance |
az.summary(trace_10_6, credible_interval=.89, round_to=2)
mean | sd | hpd_5.5% | hpd_94.5% | mcse_mean | mcse_sd | ess_mean | ess_sd | ess_bulk | ess_tail | r_hat | |
---|---|---|---|---|---|---|---|---|---|---|---|
a | -0.83 | 0.05 | -0.91 | -0.75 | 0.0 | 0.0 | 545.93 | 541.82 | 557.31 | 712.82 | 1.0 |
bm | 0.61 | 0.06 | 0.50 | 0.70 | 0.0 | 0.0 | 584.23 | 576.80 | 590.43 | 687.70 | 1.0 |
post = pm.trace_to_dataframe(trace_10_6)
p_admit_male = logistic(post['a'] + post['bm'])
p_admit_female = logistic(post['a'])
diff_admit = p_admit_male - p_admit_female
diff_admit.describe(percentiles=[.025, .5, .975])[['2.5%', '50%', '97.5%']]
2.5% 0.113468 50% 0.141427 97.5% 0.170919 dtype: float64
for i in range(6):
x = 1 + 2 * i
y1 = d_ad.admit[x] / d_ad.applications[x]
y2 = d_ad.admit[x+1] / d_ad.applications[x+1]
plt.plot([x, x+1], [y1, y2], '-C0o', lw=2)
plt.text(x + 0.25, (y1+y2)/2 + 0.05, d_ad.dept[x])
plt.ylim(0, 1);
d_ad['dept_id'] = pd.Categorical(d_ad['dept']).codes
with pm.Model() as model_10_8:
a = pm.Normal('a', 0, 10, shape=len(d_ad['dept'].unique()))
p = pm.math.invlogit(a[d_ad['dept_id'].values])
admit = pm.Binomial('admit', p=p, n=d_ad['applications'], observed=d_ad['admit'])
trace_10_8 = pm.sample(1000, tune=1000)
with pm.Model() as model_10_9:
a = pm.Normal('a', 0, 10, shape=len(d_ad['dept'].unique()))
bm = pm.Normal('bm', 0, 10)
p = pm.math.invlogit(a[d_ad['dept_id'].values] + bm * d_ad['male'])
admit = pm.Binomial('admit', p=p, n=d_ad['applications'], observed=d_ad['admit'])
trace_10_9 = pm.sample(1000, tune=1000)
Auto-assigning NUTS sampler... Initializing NUTS using jitter+adapt_diag... Multiprocess sampling (2 chains in 2 jobs) NUTS: [a]
HBox(children=(IntProgress(value=0, description='Sampling 2 chains', max=4000, style=ProgressStyle(description…
Auto-assigning NUTS sampler... Initializing NUTS using jitter+adapt_diag... Multiprocess sampling (2 chains in 2 jobs) NUTS: [bm, a]
HBox(children=(IntProgress(value=0, description='Sampling 2 chains', max=4000, style=ProgressStyle(description…
# WAIC values still off
# Plus warning flag
comp_df = az.compare({'m10.6' : trace_10_6,
'm10.7' : trace_10_7,
'm10.8' : trace_10_8,
'm10.9' : trace_10_9})
comp_df
/home/osvaldo/proyectos/00_PyMC3/arviz/arviz/stats/stats.py:1078: UserWarning: For one or more samples the posterior variance of the log predictive densities exceeds 0.4. This could be indication of WAIC starting to fail see http://arxiv.org/abs/1507.04544 for details """ /home/osvaldo/proyectos/00_PyMC3/arviz/arviz/stats/stats.py:1078: UserWarning: For one or more samples the posterior variance of the log predictive densities exceeds 0.4. This could be indication of WAIC starting to fail see http://arxiv.org/abs/1507.04544 for details """ /home/osvaldo/proyectos/00_PyMC3/arviz/arviz/stats/stats.py:1078: UserWarning: For one or more samples the posterior variance of the log predictive densities exceeds 0.4. This could be indication of WAIC starting to fail see http://arxiv.org/abs/1507.04544 for details """ /home/osvaldo/proyectos/00_PyMC3/arviz/arviz/stats/stats.py:1078: UserWarning: For one or more samples the posterior variance of the log predictive densities exceeds 0.4. This could be indication of WAIC starting to fail see http://arxiv.org/abs/1507.04544 for details """
waic | p_waic | d_waic | weight | se | dse | warning | waic_scale | |
---|---|---|---|---|---|---|---|---|
m10.8 | 104.559 | 6.25291 | 0 | 0.825254 | 304.24 | 0 | True | deviance |
m10.9 | 109.266 | 9.82074 | 4.70687 | 0.174746 | 308.669 | 3.78812 | True | deviance |
m10.6 | 990.347 | 111.824 | 885.788 | 1.14087e-64 | 16.083 | 309.771 | True | deviance |
m10.7 | 1041.33 | 79.544 | 936.767 | 4.64511e-64 | 14.6071 | 309.681 | True | deviance |
az.summary(trace_10_9, credible_interval=.89, round_to=2)
mean | sd | hpd_5.5% | hpd_94.5% | mcse_mean | mcse_sd | ess_mean | ess_sd | ess_bulk | ess_tail | r_hat | |
---|---|---|---|---|---|---|---|---|---|---|---|
a[0] | 0.68 | 0.10 | 0.52 | 0.85 | 0.0 | 0.0 | 1349.11 | 1343.82 | 1337.02 | 1459.51 | 1.0 |
a[1] | 0.64 | 0.12 | 0.46 | 0.83 | 0.0 | 0.0 | 1533.50 | 1513.36 | 1528.06 | 1194.52 | 1.0 |
a[2] | -0.58 | 0.08 | -0.70 | -0.47 | 0.0 | 0.0 | 1641.06 | 1569.64 | 1639.48 | 1492.21 | 1.0 |
a[3] | -0.61 | 0.09 | -0.74 | -0.47 | 0.0 | 0.0 | 1583.56 | 1558.25 | 1583.39 | 1528.38 | 1.0 |
a[4] | -1.06 | 0.10 | -1.23 | -0.91 | 0.0 | 0.0 | 2534.98 | 2522.44 | 2533.95 | 1494.97 | 1.0 |
a[5] | -2.63 | 0.15 | -2.89 | -2.41 | 0.0 | 0.0 | 2696.22 | 2671.08 | 2715.85 | 1530.65 | 1.0 |
bm | -0.10 | 0.08 | -0.22 | 0.03 | 0.0 | 0.0 | 1140.35 | 1128.16 | 1123.82 | 1265.51 | 1.0 |
Replicated model above but with MCMC in book.
import statsmodels.api as sm
from patsy import dmatrix
endog = d_ad.loc[:,['admit', 'reject']].values # cbind(admit,reject)
m10_7glm = sm.GLM(endog, dmatrix('~ 1', data=d_ad),
family=sm.families.Binomial())
m10_6glm = sm.GLM(endog, dmatrix('~ male', data=d_ad),
family=sm.families.Binomial())
m10_8glm = sm.GLM(endog, dmatrix('~ dept_id', data=d_ad),
family=sm.families.Binomial())
m10_9glm = sm.GLM(endog, dmatrix('~ male + dept_id', data=d_ad),
family=sm.families.Binomial())
# res = m10_7glm.fit()
# res.summary()
import statsmodels.formula.api as smf
m10_4glm = smf.glm(formula='pulled_left ~ actor + prosoc_left*condition - condition', data=d,
family=sm.families.Binomial())
pm.GLM.from_formula('pulled_left ~ actor + prosoc_left*condition - condition',
family='binomial', data=d)
# outcome and predictor almost perfectly associated
y = np.hstack([np.ones(10,)*0, np.ones(10,)])
x = np.hstack([np.ones(9,)*-1, np.ones(11,)])
m_bad = smf.glm(formula='y ~ x',
data=pd.DataFrame({'y':y, 'x':x}),
family=sm.families.Binomial()).fit()
m_bad.summary()
Dep. Variable: | y | No. Observations: | 20 |
---|---|---|---|
Model: | GLM | Df Residuals: | 18 |
Model Family: | Binomial | Df Model: | 1 |
Link Function: | logit | Scale: | 1.0000 |
Method: | IRLS | Log-Likelihood: | -3.3510 |
Date: | Sat, 29 Jun 2019 | Deviance: | 6.7020 |
Time: | 12:49:51 | Pearson chi2: | 11.0 |
No. Iterations: | 21 | Covariance Type: | nonrobust |
coef | std err | z | P>|z| | [0.025 | 0.975] | |
---|---|---|---|---|---|---|
Intercept | -10.1317 | 8032.690 | -0.001 | 0.999 | -1.58e+04 | 1.57e+04 |
x | 12.4343 | 8032.690 | 0.002 | 0.999 | -1.57e+04 | 1.58e+04 |
with pm.Model() as m_good:
ab = pm.Normal('ab', 0, 10, shape=2)
p = pm.math.invlogit(ab[0] + ab[1]*x)
y_ = pm.Binomial('y_', 1, p, observed=y)
MAP = pm.find_MAP()
MAP
/home/osvaldo/proyectos/00_PyMC3/pymc3/pymc3/tuning/starting.py:61: UserWarning: find_MAP should not be used to initialize the NUTS sampler, simply call pymc3.sample() and it will automatically initialize NUTS in a better way. warnings.warn('find_MAP should not be used to initialize the NUTS sampler, simply call pymc3.sample() and it will automatically initialize NUTS in a better way.')
HBox(children=(IntProgress(value=0, max=5000), HTML(value='')))
{'ab': array([-1.72704484, 4.01710522])}
trace = pm.sample(1000, tune=1000, model=m_good)
tracedf = pm.trace_to_dataframe(trace)
grid = (sns.PairGrid(tracedf,
diag_sharey=False)
.map_diag(sns.kdeplot)
.map_upper(plt.scatter, alpha=0.1))
Auto-assigning NUTS sampler... Initializing NUTS using jitter+adapt_diag... Multiprocess sampling (2 chains in 2 jobs) NUTS: [ab]
HBox(children=(IntProgress(value=0, description='Sampling 2 chains', max=4000, style=ProgressStyle(description…
There were 141 divergences after tuning. Increase `target_accept` or reparameterize. The acceptance probability does not match the target. It is 0.45368539279407266, but should be close to 0.8. Try to increase the number of tuning steps. The acceptance probability does not match the target. It is 0.8787129508505335, but should be close to 0.8. Try to increase the number of tuning steps. The estimated number of effective samples is smaller than 200 for some parameters.
y = stats.binom.rvs(n=1000, p=1/1000, size=100000)
np.mean(y), np.var(y)
(0.99678, 1.0015896316)
dk = pd.read_csv('Data/Kline', sep=';')
dk
culture | population | contact | total_tools | mean_TU | |
---|---|---|---|---|---|
0 | Malekula | 1100 | low | 13 | 3.2 |
1 | Tikopia | 1500 | low | 22 | 4.7 |
2 | Santa Cruz | 3600 | low | 24 | 4.0 |
3 | Yap | 4791 | high | 43 | 5.0 |
4 | Lau Fiji | 7400 | high | 33 | 5.0 |
5 | Trobriand | 8000 | high | 19 | 4.0 |
6 | Chuuk | 9200 | high | 40 | 3.8 |
7 | Manus | 13000 | low | 28 | 6.6 |
8 | Tonga | 17500 | high | 55 | 5.4 |
9 | Hawaii | 275000 | low | 71 | 6.6 |
dk.log_pop = np.log(dk.population)
dk.contact_high = (dk.contact == "high").astype(int)
/home/osvaldo/anaconda3/lib/python3.7/site-packages/ipykernel_launcher.py:1: UserWarning: Pandas doesn't allow columns to be created via a new attribute name - see https://pandas.pydata.org/pandas-docs/stable/indexing.html#attribute-access """Entry point for launching an IPython kernel. /home/osvaldo/anaconda3/lib/python3.7/site-packages/ipykernel_launcher.py:2: UserWarning: Pandas doesn't allow columns to be created via a new attribute name - see https://pandas.pydata.org/pandas-docs/stable/indexing.html#attribute-access
from theano import shared
# casting data to theano shared variable.
# It is for out of sample prediction from model with sampled trace
log_pop = shared(dk.log_pop.values)
contact_high = shared(dk.contact_high.values)
total_tools = shared(dk.total_tools.values)
with pm.Model() as m_10_10:
a = pm.Normal('a', 0, 100)
b = pm.Normal('b', 0, 1, shape=3)
lam = pm.math.exp(a + b[0] * log_pop + b[1] * contact_high + b[2] * contact_high * log_pop)
obs = pm.Poisson('total_tools', lam, observed=total_tools)
trace_10_10 = pm.sample(1000, tune=1000)
Auto-assigning NUTS sampler... Initializing NUTS using jitter+adapt_diag... Multiprocess sampling (2 chains in 2 jobs) NUTS: [b, a]
HBox(children=(IntProgress(value=0, description='Sampling 2 chains', max=4000, style=ProgressStyle(description…
summary = az.summary(trace_10_10, credible_interval=.89)[['mean', 'sd', 'hpd_5.5%', 'hpd_94.5%']]
trace_cov = pm.trace_cov(trace_10_10, model=m_10_10)
invD = (np.sqrt(np.diag(trace_cov))**-1)[:, None]
trace_corr = pd.DataFrame(invD*trace_cov*invD.T, index=summary.index, columns=summary.index)
summary.join(trace_corr).round(2)
mean | sd | hpd_5.5% | hpd_94.5% | a | b[0] | b[1] | b[2] | |
---|---|---|---|---|---|---|---|---|
a | 0.93 | 0.36 | 0.34 | 1.49 | 1.00 | -0.98 | -0.03 | -0.03 |
b[0] | 0.26 | 0.04 | 0.21 | 0.32 | -0.98 | 1.00 | 0.03 | 0.01 |
b[1] | -0.09 | 0.84 | -1.48 | 1.20 | -0.03 | 0.03 | 1.00 | -0.99 |
b[2] | 0.04 | 0.09 | -0.11 | 0.19 | -0.03 | 0.01 | -0.99 | 1.00 |
az.plot_forest(trace_10_10);
lambda_high = np.exp(trace_10_10['a'] + trace_10_10['b'][:,1] + (trace_10_10['b'][:,0] + trace_10_10['b'][:,2]) * 8)
lambda_low = np.exp(trace_10_10['a'] + trace_10_10['b'][:,0] * 8 )
diff = lambda_high - lambda_low
np.sum(diff > 0) / len(diff)
0.9625
with pm.Model() as m_10_11:
a = pm.Normal('a', 0, 100)
b = pm.Normal('b', 0, 1, shape=2)
lam = pm.math.exp(a + b[0] * log_pop + b[1] * contact_high)
obs = pm.Poisson('total_tools', lam, observed=total_tools)
trace_10_11 = pm.sample(1000, tune=1000)
Auto-assigning NUTS sampler... Initializing NUTS using jitter+adapt_diag... Multiprocess sampling (2 chains in 2 jobs) NUTS: [b, a]
HBox(children=(IntProgress(value=0, description='Sampling 2 chains', max=4000, style=ProgressStyle(description…
with pm.Model() as m_10_12:
a = pm.Normal('a', 0, 100)
b = pm.Normal('b', 0, 1)
lam = pm.math.exp(a + b * log_pop)
obs = pm.Poisson('total_tools', lam, observed=total_tools)
trace_10_12 = pm.sample(1000, tune=1000)
with pm.Model() as m_10_13:
a = pm.Normal('a', 0, 100)
b = pm.Normal('b', 0, 1)
lam = pm.math.exp(a + b * contact_high)
obs = pm.Poisson('total_tools', lam, observed=total_tools)
trace_10_13 = pm.sample(1000, tune=1000)
Auto-assigning NUTS sampler... Initializing NUTS using jitter+adapt_diag... Multiprocess sampling (2 chains in 2 jobs) NUTS: [b, a]
HBox(children=(IntProgress(value=0, description='Sampling 2 chains', max=4000, style=ProgressStyle(description…
The acceptance probability does not match the target. It is 0.8835601127234468, but should be close to 0.8. Try to increase the number of tuning steps. The number of effective samples is smaller than 25% for some parameters.
Auto-assigning NUTS sampler... Initializing NUTS using jitter+adapt_diag... Multiprocess sampling (2 chains in 2 jobs) NUTS: [b, a]
HBox(children=(IntProgress(value=0, description='Sampling 2 chains', max=4000, style=ProgressStyle(description…
with pm.Model() as m_10_14:
a = pm.Normal('a', 0, 100)
lam = pm.math.exp(a)
obs = pm.Poisson('total_tools', lam, observed=total_tools)
trace_10_14 = pm.sample(1000, tune=1000)
Auto-assigning NUTS sampler... Initializing NUTS using jitter+adapt_diag... Multiprocess sampling (2 chains in 2 jobs) NUTS: [a]
HBox(children=(IntProgress(value=0, description='Sampling 2 chains', max=4000, style=ProgressStyle(description…
traces = [trace_10_10, trace_10_11, trace_10_12, trace_10_13, trace_10_14]
models = [m_10_10, m_10_11, m_10_12, m_10_13, m_10_14]
model_names = ['m10.10', 'm10.11', 'm10.12', 'm10.13','m10.14']
dictionary = dict(zip(model_names, traces))
islands_compare = az.compare(dictionary)
islands_compare
/home/osvaldo/proyectos/00_PyMC3/arviz/arviz/stats/stats.py:1078: UserWarning: For one or more samples the posterior variance of the log predictive densities exceeds 0.4. This could be indication of WAIC starting to fail see http://arxiv.org/abs/1507.04544 for details """ /home/osvaldo/proyectos/00_PyMC3/arviz/arviz/stats/stats.py:1078: UserWarning: For one or more samples the posterior variance of the log predictive densities exceeds 0.4. This could be indication of WAIC starting to fail see http://arxiv.org/abs/1507.04544 for details """ /home/osvaldo/proyectos/00_PyMC3/arviz/arviz/stats/stats.py:1078: UserWarning: For one or more samples the posterior variance of the log predictive densities exceeds 0.4. This could be indication of WAIC starting to fail see http://arxiv.org/abs/1507.04544 for details """ /home/osvaldo/proyectos/00_PyMC3/arviz/arviz/stats/stats.py:1078: UserWarning: For one or more samples the posterior variance of the log predictive densities exceeds 0.4. This could be indication of WAIC starting to fail see http://arxiv.org/abs/1507.04544 for details """ /home/osvaldo/proyectos/00_PyMC3/arviz/arviz/stats/stats.py:1078: UserWarning: For one or more samples the posterior variance of the log predictive densities exceeds 0.4. This could be indication of WAIC starting to fail see http://arxiv.org/abs/1507.04544 for details """
waic | p_waic | d_waic | weight | se | dse | warning | waic_scale | |
---|---|---|---|---|---|---|---|---|
m10.11 | 79.355 | 4.39772 | 0 | 0.503126 | 11.2484 | 0 | True | deviance |
m10.10 | 80.6378 | 5.15511 | 1.28279 | 0.276046 | 10.9265 | 1.18494 | True | deviance |
m10.12 | 84.4378 | 3.74629 | 5.08275 | 0.219833 | 8.54095 | 7.74886 | True | deviance |
m10.14 | 141.562 | 8.28172 | 62.2074 | 0.000230229 | 43.6218 | 32.9269 | True | deviance |
m10.13 | 150.564 | 16.7978 | 71.2093 | 0.000765188 | 30.4152 | 44.7903 | True | deviance |
az.plot_compare(islands_compare);
# set new value for out-of-sample prediction
log_pop_seq = np.linspace(6, 13, 30)
log_pop.set_value(np.hstack([log_pop_seq, log_pop_seq]))
contact_high.set_value(np.hstack([np.repeat(0, 30), np.repeat(1, 30)]))
islands_ensemble = pm.sample_posterior_predictive_w(traces, 10000, models,
weights=islands_compare.weight.sort_index(ascending=True))
HBox(children=(IntProgress(value=0, max=10000), HTML(value='')))
_, axes = plt.subplots(1, 1, figsize=(5, 5))
index = dk.contact_high==1
axes.scatter(np.log(dk.population)[~index], dk.total_tools[~index],
facecolors='none', edgecolors='k', lw=1)
axes.scatter(np.log(dk.population)[index], dk.total_tools[index])
mp = islands_ensemble['total_tools'][:, :30]
axes.plot(log_pop_seq, np.median(mp, axis=0), '--', color='k')
az.plot_hpd(log_pop_seq, mp, credible_interval=.5, color='k')
mp = islands_ensemble['total_tools'][:, 30:]
axes.plot(log_pop_seq, np.median(mp, axis=0))
az.plot_hpd(log_pop_seq, mp, credible_interval=.5)
axes.set_xlabel('log-population')
axes.set_ylabel('total tools')
axes.set_xlim(6.8, 12.8)
axes.set_ylim(10, 73);
This is the same as 10.41, but in the book using MCMC rather than MAP.
log_pop_c = dk.log_pop.values - dk.log_pop.values.mean()
log_pop.set_value(log_pop_c)
contact_high.set_value(dk.contact_high.values)
total_tools.set_value(dk.total_tools.values)
with pm.Model() as m_10_10c:
a = pm.Normal('a', 0, 100)
b = pm.Normal('b', 0, 1, shape=3)
lam = pm.math.exp(a + b[0] * log_pop + b[1] * contact_high + b[2] * contact_high * log_pop)
obs = pm.Poisson('total_tools', lam, observed=total_tools)
trace_10_10c = pm.sample(1000, tune=1000)
Auto-assigning NUTS sampler... Initializing NUTS using jitter+adapt_diag... Multiprocess sampling (2 chains in 2 jobs) NUTS: [b, a]
HBox(children=(IntProgress(value=0, description='Sampling 2 chains', max=4000, style=ProgressStyle(description…
az.summary(trace_10_10c, credible_interval=.89, round_to=2)
mean | sd | hpd_5.5% | hpd_94.5% | mcse_mean | mcse_sd | ess_mean | ess_sd | ess_bulk | ess_tail | r_hat | |
---|---|---|---|---|---|---|---|---|---|---|---|
a | 3.31 | 0.09 | 3.19 | 3.46 | 0.0 | 0.0 | 1218.90 | 1216.39 | 1217.22 | 1218.70 | 1.0 |
b[0] | 0.26 | 0.04 | 0.21 | 0.33 | 0.0 | 0.0 | 1571.34 | 1560.62 | 1578.04 | 1140.98 | 1.0 |
b[1] | 0.28 | 0.12 | 0.10 | 0.47 | 0.0 | 0.0 | 1292.44 | 1277.52 | 1293.50 | 1549.39 | 1.0 |
b[2] | 0.06 | 0.17 | -0.21 | 0.32 | 0.0 | 0.0 | 1409.02 | 1272.78 | 1408.09 | 1396.35 | 1.0 |
for trace in [trace_10_10, trace_10_10c]:
tracedf = pm.trace_to_dataframe(trace)
grid = (sns.PairGrid(tracedf,
diag_sharey=False)
.map_diag(sns.kdeplot)
.map_upper(plt.scatter, alpha=0.1))
num_days = 30
y = np.random.poisson(1.5, num_days)
num_weeks = 4
y_new = np.random.poisson(0.5*7, num_weeks)
y_all = np.hstack([y, y_new])
exposure = np.hstack([np.repeat(1, 30), np.repeat(7, 4)]).astype('float')
monastery = np.hstack([np.repeat(0, 30), np.repeat(1, 4)])
log_days = np.log(exposure)
with pm.Model() as m_10_15:
a = pm.Normal('a', 0., 100.)
b = pm.Normal('b', 0., 1.)
lam = pm.math.exp(log_days + a + b*monastery)
obs = pm.Poisson('y', lam, observed=y_all)
trace_10_15 = pm.sample(1000, tune=1000)
Auto-assigning NUTS sampler... Initializing NUTS using jitter+adapt_diag... Multiprocess sampling (2 chains in 2 jobs) NUTS: [b, a]
HBox(children=(IntProgress(value=0, description='Sampling 2 chains', max=4000, style=ProgressStyle(description…
trace_10_15.add_values(dict(lambda_old=np.exp(trace_10_15['a']),
lambda_new=np.exp(trace_10_15['a'] + trace_10_15['b'])))
az.summary(trace_10_15, var_names=['lambda_old', 'lambda_new'],
credible_interval=.89, round_to=2)
mean | sd | hpd_5.5% | hpd_94.5% | mcse_mean | mcse_sd | ess_mean | ess_sd | ess_bulk | ess_tail | r_hat | |
---|---|---|---|---|---|---|---|---|---|---|---|
lambda_old | 1.95 | 0.25 | 1.54 | 2.34 | 0.01 | 0.0 | 1474.08 | 1474.08 | 1455.16 | 1299.94 | 1.0 |
lambda_new | 0.51 | 0.13 | 0.29 | 0.70 | 0.00 | 0.0 | 1589.43 | 1589.43 | 1525.23 | 1421.27 | 1.0 |
# simulate career choices among 500 individuals
N = 500 # number of individuals
income = np.arange(3) + 1 # expected income of each career
score = 0.5 * income # scores for each career, based on income
# next line converts scores to probabilities
def softmax(w):
e = np.exp(w)
return e/np.sum(e, axis=0)
p = softmax(score)
# now simulate choice
# outcome career holds event type values, not counts
career = np.random.multinomial(1, p, size=N)
career = np.where(career==1)[1]
career[:11]
array([1, 0, 1, 2, 0, 1, 0, 1, 2, 2, 2])
with pm.Model() as m_10_16:
b = pm.Normal('b', 0., 5.)
s2 = b*2
s3 = b*3
p_ = pm.math.stack([0, s2, s3])
obs = pm.Categorical('career', p=tt.nnet.softmax(p_), observed=career)
trace_10_16 = pm.sample(1000, tune=2000, cores=2)
az.summary(trace_10_16, credible_interval=.89, round_to=2)
/anaconda/envs/fund/lib/python3.7/site-packages/ipykernel_launcher.py:7: UserWarning: DEPRECATION: If x is a vector, Softmax will not automatically pad x anymore in next releases. If you need it, please do it manually. The vector case is gonna be supported soon and the output will be a vector. import sys Auto-assigning NUTS sampler... Initializing NUTS using jitter+adapt_diag... Multiprocess sampling (2 chains in 2 jobs) NUTS: [b] Sampling 2 chains, 0 divergences: 100%|██████████| 6000/6000 [00:11<00:00, 530.57draws/s]
mean | sd | hpd_5.5% | hpd_94.5% | mcse_mean | mcse_sd | ess_mean | ess_sd | ess_bulk | ess_tail | r_hat | |
---|---|---|---|---|---|---|---|---|---|---|---|
b | 0.33 | 0.04 | 0.27 | 0.4 | 0.0 | 0.0 | 615.94 | 615.94 | 623.48 | 1125.77 | 1.01 |
N = 100
# simulate family incomes for each individual
family_income = np.random.rand(N)
# assign a unique coefficient for each type of event
b = np.arange(3)-1
p = softmax(score[:, None] + np.outer(b, family_income)).T
career = np.asarray([np.random.multinomial(1, pp) for pp in p])
career = np.where(career==1)[1]
career
array([1, 2, 2, 2, 1, 2, 2, 2, 0, 2, 2, 2, 1, 2, 2, 2, 1, 2, 1, 2, 0, 2, 2, 0, 2, 1, 2, 2, 2, 2, 1, 2, 2, 1, 2, 0, 2, 1, 2, 2, 2, 1, 2, 1, 1, 2, 1, 2, 2, 1, 1, 1, 2, 2, 2, 1, 1, 2, 2, 2, 2, 1, 2, 0, 1, 2, 2, 2, 2, 2, 2, 1, 0, 1, 1, 0, 2, 2, 0, 2, 2, 2, 2, 0, 2, 2, 2, 2, 2, 1, 2, 1, 2, 2, 2, 2, 1, 2, 2, 0])
with pm.Model() as m_10_17:
a23 = pm.Normal('a23', 0., 5., shape=2)
b23 = pm.Normal('b23', 0., 5., shape=2)
s2 = a23[0] + b23[0]*family_income
s3 = a23[1] + b23[1]*family_income
p_ = pm.math.stack([np.zeros(N), s2, s3]).T
obs = pm.Categorical('career', p=tt.nnet.softmax(p_), observed=career)
trace_10_17 = pm.sample(1000, tune=2000, cores=2)
az.summary(trace_10_17, credible_interval=.89, round_to=2)
Auto-assigning NUTS sampler... Initializing NUTS using jitter+adapt_diag... Multiprocess sampling (2 chains in 2 jobs) NUTS: [b23, a23] Sampling 2 chains, 0 divergences: 100%|██████████| 6000/6000 [00:38<00:00, 157.53draws/s]
mean | sd | hpd_5.5% | hpd_94.5% | mcse_mean | mcse_sd | ess_mean | ess_sd | ess_bulk | ess_tail | r_hat | |
---|---|---|---|---|---|---|---|---|---|---|---|
a23[0] | 0.87 | 0.76 | -0.37 | 2.03 | 0.03 | 0.02 | 555.77 | 555.77 | 557.28 | 788.43 | 1.0 |
a23[1] | 1.76 | 0.68 | 0.70 | 2.86 | 0.03 | 0.02 | 549.78 | 527.32 | 557.30 | 691.57 | 1.0 |
b23[0] | 0.27 | 1.27 | -1.99 | 2.07 | 0.05 | 0.04 | 569.73 | 569.73 | 570.80 | 842.66 | 1.0 |
b23[1] | 0.34 | 1.15 | -1.32 | 2.29 | 0.05 | 0.03 | 560.08 | 560.08 | 557.17 | 834.00 | 1.0 |
d_ad = pd.read_csv('./Data/UCBadmit.csv', sep=';')
# binomial model of overall admission probability
with pm.Model() as m_binom:
a = pm.Normal('a', 0, 100)
p = pm.math.invlogit(a)
admit = pm.Binomial('admit', p=p, n=d_ad.applications, observed=d_ad.admit)
trace_binom = pm.sample(1000, tune=1000)
# Poisson model of overall admission rate and rejection rate
with pm.Model() as m_pois:
a = pm.Normal('a', 0, 100, shape=2)
lam = pm.math.exp(a)
admit = pm.Poisson('admit', lam[0], observed=d_ad.admit)
rej = pm.Poisson('rej', lam[1], observed=d_ad.reject)
trace_pois = pm.sample(1000, tune=1000)
Auto-assigning NUTS sampler... Initializing NUTS using jitter+adapt_diag... Multiprocess sampling (2 chains in 2 jobs) NUTS: [a]
HBox(children=(IntProgress(value=0, description='Sampling 2 chains', max=4000, style=ProgressStyle(description…
Auto-assigning NUTS sampler... Initializing NUTS using jitter+adapt_diag... Multiprocess sampling (2 chains in 2 jobs) NUTS: [a]
HBox(children=(IntProgress(value=0, description='Sampling 2 chains', max=4000, style=ProgressStyle(description…
m_binom = pm.summary(trace_binom, alpha=.11).round(2)
logistic(m_binom['mean'])
a 0.386986 Name: mean, dtype: float64
m_pois = pm.summary(trace_pois, alpha=.11).round(2)
m_pois['mean'][0]
np.exp(m_pois['mean'][0])/(np.exp(m_pois['mean'][0])+np.exp(m_pois['mean'][1]))
0.38936076605077796
# simulate
N = 100
x = np.random.rand(N)
y = np.random.geometric(logistic(-1 + 2*x), size=N)
with pm.Model() as m_10_18:
a = pm.Normal('a', 0, 10)
b = pm.Normal('b', 0, 1)
p = pm.math.invlogit(a + b*x)
obs = pm.Geometric('y', p=p, observed=y)
trace_10_18 = pm.sample(1000, tune=1000)
az.summary(trace_10_18, credible_interval=.89, round_to=2)
Auto-assigning NUTS sampler... Initializing NUTS using jitter+adapt_diag... Multiprocess sampling (2 chains in 2 jobs) NUTS: [b, a]
HBox(children=(IntProgress(value=0, description='Sampling 2 chains', max=4000, style=ProgressStyle(description…
mean | sd | hpd_5.5% | hpd_94.5% | mcse_mean | mcse_sd | ess_mean | ess_sd | ess_bulk | ess_tail | r_hat | |
---|---|---|---|---|---|---|---|---|---|---|---|
a | -0.98 | 0.22 | -1.33 | -0.64 | 0.01 | 0.01 | 644.24 | 635.21 | 646.45 | 666.90 | 1.0 |
b | 1.85 | 0.44 | 1.16 | 2.55 | 0.02 | 0.01 | 662.98 | 662.52 | 664.11 | 817.05 | 1.0 |
import sys, IPython, scipy, matplotlib, platform
print("""This notebook was created using:\nPython {}\nIPython {}\nPyMC3 {}\nArviZ {}\nNumPy {}\nSciPy {}\nMatplotlib {}\n""".format(sys.version[:5], IPython.__version__, pm.__version__, az.__version__, np.__version__, scipy.__version__, matplotlib.__version__))
This notebook was created using: Python 3.7.3 IPython 7.5.0 PyMC3 3.7 ArviZ 0.4.1 NumPy 1.16.4 SciPy 1.2.1 Matplotlib 3.1.0