import numpy as np, pandas as pd, pymc as pm, pytensor.tensor as pt, arviz as az
from scipy.special import expit
from seeds_gibbs import seeds_gibbs
d = pd.read_csv('surgical.csv'); r = d['r'].to_numpy(); n = d['n'].to_numpy(); I = len(d)
with pm.Model() as m:
    mu = pm.Normal('mu', 0, 10); sigma = pm.HalfNormal('sigma', 1.0)
    z = pm.Normal('z', 0, 1, shape=I); b = pm.Deterministic('b', sigma*z)
    p = pm.Deterministic('p', pm.math.invlogit(mu + b))
    pm.Binomial('r', n=n, p=p, observed=r)
    idata = pm.sample(2000, tune=2000, chains=4, target_accept=0.95, random_seed=1, progressbar=False)
po = idata.posterior; P_pm = po['p'].values.reshape(-1, I)
print('PyMC: population rate %.1f%%  sigma %.3f  max r_hat %.3f'
      % (100*expit(po['mu'].values).mean(), float(po['sigma'].mean()), float(az.summary(idata,var_names=['mu','sigma'])['r_hat'].max())))
g = seeds_gibbs(r.astype(float), n.astype(float), np.ones((I,1)), R=30000, burn=8000, seed=1)
P_fs = expit(g['beta'][:,0][:,None] + g['bdraws'])
print('\nhosp  raw%   PyMC%   from-scratch%')
for i in range(I): print('  %2d  %5.1f  %6.1f  %12.1f' % (i+1, 100*r[i]/n[i], 100*P_pm[:,i].mean(), 100*P_fs[:,i].mean()))

Initializing NUTS using jitter+adapt_diag...

Multiprocess sampling (4 chains in 4 jobs)

NUTS: [mu, sigma, z]

Sampling 4 chains for 2_000 tune and 2_000 draw iterations (8_000 + 8_000 draws total) took 5 seconds.

PyMC: population rate 7.2%  sigma 0.445  max r_hat 1.000

hosp  raw%   PyMC%   from-scratch%
   1    0.0     5.1           5.8
   2   12.2    10.6           9.9
   3    6.7     7.0           7.2
   4    5.7     5.9           6.1
   5    3.8     5.0           5.6
   6    6.6     6.9           7.1
   7    6.1     6.6           6.9
   8   14.4    12.6          11.6
   9    6.8     7.0           7.1
  10    8.2     7.9           7.9
  11   11.3    10.4           9.8
  12    6.7     6.8           7.0

Surgical — PyMC cross-check (intercept-only binomial GLMM)¶

Results¶