import numpy as np, pandas as pd, pymc as pm, pytensor.tensor as pt, arviz as az
from scipy.special import expit
from betabin_gibbs import betabin_gibbs
d = pd.read_csv('baseball.csv'); r = d['r'].to_numpy(); n = d['n'].to_numpy(); I = len(d); truth = d['remaining_avg'].to_numpy()
# beta-binomial
with pm.Model() as mbb:
    mu = pm.Beta('mu', 1, 1); kappa = pm.HalfCauchy('kappa', 50)
    a = mu*kappa; b = (1-mu)*kappa
    pm.BetaBinomial('r', alpha=a, beta=b, n=n, observed=r)
    pm.Deterministic('p', (a + r) / (a + b + n))               # shrunk posterior-mean rates
    idb = pm.sample(2000, tune=2000, chains=4, target_accept=0.95, random_seed=1, progressbar=False)
# logit-normal
with pm.Model() as mln:
    m0 = pm.Normal('m0', 0, 5); sg = pm.HalfNormal('sg', 1.0); z = pm.Normal('z', 0, 1, shape=I)
    p = pm.Deterministic('p', pm.math.invlogit(m0 + sg*z))
    pm.Binomial('r', n=n, p=p, observed=r)
    idl = pm.sample(2000, tune=2000, chains=4, target_accept=0.95, random_seed=1, progressbar=False)
pbb = idb.posterior['p'].values.reshape(-1,I).mean(0); pln = idl.posterior['p'].values.reshape(-1,I).mean(0)
g = betabin_gibbs(r.astype(float), n.astype(float), R=30000, burn=8000, seed=1)['p'].mean(0)
def rmse(x): return np.sqrt(np.mean((x-truth)**2))
print('PyMC beta-binomial: mu %.3f kappa(med) %.0f   RMSE %.4f'%(float(idb.posterior['mu'].mean()),float(np.median(idb.posterior['kappa'].values)),rmse(pbb)))
print('PyMC logit-normal : sigma %.3f               RMSE %.4f'%(float(idl.posterior['sg'].mean()),rmse(pln)))
print('raw RMSE %.4f | beta-binom vs logit-normal corr %.4f | PyMC-vs-fromscratch(betabin) max diff %.4f'
      %(rmse(r/n), np.corrcoef(pbb,pln)[0,1], np.abs(pbb-g).max()))

Initializing NUTS using jitter+adapt_diag...

Multiprocess sampling (4 chains in 4 jobs)

NUTS: [mu, kappa]

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

C:\Users\user\anaconda3\envs\pymc-env\Lib\site-packages\threadpoolctl.py:1226: RuntimeWarning: 
Found Intel OpenMP ('libiomp') and LLVM OpenMP ('libomp') loaded at
the same time. Both libraries are known to be incompatible and this
can cause random crashes or deadlocks on Linux when loaded in the
same Python program.
Using threadpoolctl may cause crashes or deadlocks. For more
information and possible workarounds, please see
    https://github.com/joblib/threadpoolctl/blob/master/multiple_openmp.md

  warnings.warn(msg, RuntimeWarning)
Initializing NUTS using jitter+adapt_diag...

Multiprocess sampling (4 chains in 4 jobs)

NUTS: [m0, sg, z]

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

PyMC beta-binomial: mu 0.266 kappa(med) 132   RMSE 0.0390
PyMC logit-normal : sigma 0.164               RMSE 0.0384
raw RMSE 0.0690 | beta-binom vs logit-normal corr 0.9996 | PyMC-vs-fromscratch(betabin) max diff 0.0095

Baseball — PyMC cross-check (beta-binomial + logit-normal)¶

Results¶