import numpy as np, pymc as pm, pytensor.tensor as pt, arviz as az
from longlogit_gibbs import toenail_data, longlogit_gibbs
y, X, g = toenail_data(); J = int(g.max()+1); k = X.shape[1]
with pm.Model() as m:
    beta = pm.Normal('beta', 0, 10, shape=k); sigma = pm.HalfNormal('sigma', 5.0)
    z = pm.Normal('z', 0, 1, shape=J); b = pm.Deterministic('b', sigma*z)
    pm.Bernoulli('y', logit_p=pt.dot(X, beta) + b[g], observed=y)
    idata = pm.sample(1500, tune=1500, chains=4, target_accept=0.95, random_seed=1, progressbar=False)
po = idata.posterior; bb = po['beta'].values.reshape(-1,k); nm = ['intercept','treat','month','treat×month']
for j in range(k): print('%-13s %6.3f [%.2f, %.2f]' % (nm[j], bb[:,j].mean(), np.percentile(bb[:,j],2.5), np.percentile(bb[:,j],97.5)))
print('sigma = %.2f [%.1f, %.1f]   max r_hat %.3f' % (float(po['sigma'].mean()), float(po['sigma'].quantile(.025)), float(po['sigma'].quantile(.975)), float(az.summary(idata,var_names=['beta','sigma'])['r_hat'].max())))
gg = longlogit_gibbs(y, X, g, R=12000, burn=3000, seed=1)
print('\nfrom-scratch PG-Gibbs:  beta %s  sigma %.2f' % (np.round(gg['beta'].mean(0),3), gg['sigma'].mean()))

Initializing NUTS using jitter+adapt_diag...

Multiprocess sampling (4 chains in 4 jobs)

NUTS: [beta, sigma, z]

Sampling 4 chains for 1_500 tune and 1_500 draw iterations (6_000 + 6_000 draws total) took 13 seconds.

intercept     -1.663 [-2.58, -0.81]
treat         -0.183 [-1.38, 0.99]
month         -0.397 [-0.49, -0.31]
treat×month   -0.141 [-0.28, -0.01]
sigma = 4.14 [3.4, 5.0]   max r_hat 1.010

from-scratch PG-Gibbs:  beta [-1.651 -0.172 -0.396 -0.141]  sigma 4.10

Toenail — PyMC cross-check (mixed logistic, non-centred)¶

Results¶