Rugby — PyMC cross-check (crossed-effects hierarchical Poisson)¶

Companion to rugby_python.ipynb (from-scratch Metropolis-within-Gibbs). Same model, sampled jointly with NUTS, with the sum-to-zero constraint written explicitly via centred deterministics. We check that the team attack/defence estimates land on top of the from-scratch sampler.

Model¶

$$\log\theta_{\text{home}}=\mu+\text{home}+\text{att}_h-\text{def}_a,\qquad \log\theta_{\text{away}}=\mu+\text{att}_a-\text{def}_h,$$ with $\text{att}_t,\text{def}_t$ drawn from $N(0,\sigma^2_\bullet)$ and re-centred to sum to zero (atts = atts_star - mean). $\sigma_{\text{att}},\sigma_{\text{def}}\sim\text{HalfNormal}(2)$; $\mu,\text{home}\sim N(0,10)$.

In [1]:
import numpy as np, pandas as pd, pymc as pm, pytensor.tensor as pt, arviz as az, matplotlib.pyplot as plt
from rugby_gibbs import rugby_data, rugby_gibbs

d = pd.read_csv('rugby.csv', index_col=0); d = d[d.year == 2014].reset_index(drop=True)
teams = sorted(set(d.home_team) | set(d.away_team)); idx = {t: i for i, t in enumerate(teams)}; T = len(teams)
hi = d.home_team.map(idx).to_numpy(); ai = d.away_team.map(idx).to_numpy()
hs = d.home_score.to_numpy(); as_ = d.away_score.to_numpy()

with pm.Model() as m:
    sd_att = pm.HalfNormal('sd_att', 2.0); sd_def = pm.HalfNormal('sd_def', 2.0)
    mu = pm.Normal('mu', 0, 10); home = pm.Normal('home', 0, 10)
    atts_s = pm.Normal('atts_s', 0, sd_att, shape=T); defs_s = pm.Normal('defs_s', 0, sd_def, shape=T)
    atts = pm.Deterministic('atts', atts_s - atts_s.mean())
    defs = pm.Deterministic('defs', defs_s - defs_s.mean())
    pm.Poisson('home_pts', pt.exp(mu + home + atts[hi] - defs[ai]), observed=hs)
    pm.Poisson('away_pts', pt.exp(mu + atts[ai] - defs[hi]), observed=as_)
    idata = pm.sample(2000, tune=2000, chains=4, target_accept=0.9, random_seed=1, progressbar=False)

po = idata.posterior; rh = az.rhat(idata)
print(f"mu   = {float(po['mu'].mean()):.3f}   r_hat {float(rh['mu']):.3f}")
print(f"home = {float(po['home'].mean()):.3f}  95% CI [{float(po['home'].quantile(.025)):.2f}, {float(po['home'].quantile(.975)):.2f}]  -> x{np.exp(float(po['home'].mean())):.2f}   r_hat {float(rh['home']):.3f}")
print(f"sd_att = {float(po['sd_att'].mean()):.3f}   sd_def = {float(po['sd_def'].mean()):.3f}")
att = po['atts'].mean(('chain','draw')).values; deff = po['defs'].mean(('chain','draw')).values
print(f"\n{'team':<10}{'attack':>9}{'defence':>9}")
for t in np.argsort(-att): print(f"{teams[t]:<10}{att[t]:>9.3f}{deff[t]:>9.3f}")

Initializing NUTS using jitter+adapt_diag...

Multiprocess sampling (4 chains in 4 jobs)

NUTS: [sd_att, sd_def, mu, home, atts_s, defs_s]

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

mu   = 2.662   r_hat 1.000
home = 0.370  95% CI [0.19, 0.55]  -> x1.45   r_hat 1.001
sd_att = 0.467   sd_def = 0.554

team         attack  defence
England       0.404    0.249
Ireland       0.181    0.557
Wales         0.129    0.141
France        0.052   -0.138
Italy        -0.235   -0.519
Scotland     -0.530   -0.291
In [2]:
# agreement with the from-scratch Gibbs
g = rugby_gibbs(*rugby_data('rugby.csv', 2014)[:4], T, R=30000, burn=5000, seed=1)
ga = g['att'].mean(0); gd = g['deff'].mean(0)
fig, ax = plt.subplots(1, 2, figsize=(11, 4.4))
ax[0].plot([-.7,.7],[-.7,.7],'k:',lw=1); ax[0].scatter(ga, att, color='steelblue', zorder=3)
for t in range(T): ax[0].annotate(teams[t],(ga[t],att[t]),fontsize=8,xytext=(4,3),textcoords='offset points')
ax[0].set_xlabel('from-scratch attack'); ax[0].set_ylabel('PyMC attack'); ax[0].set_title('Attack (on y=x)')
ax[1].plot([-.7,.7],[-.7,.7],'k:',lw=1); ax[1].scatter(gd, deff, color='seagreen', zorder=3)
for t in range(T): ax[1].annotate(teams[t],(gd[t],deff[t]),fontsize=8,xytext=(4,3),textcoords='offset points')
ax[1].set_xlabel('from-scratch defence'); ax[1].set_ylabel('PyMC defence'); ax[1].set_title('Defence (on y=x)')
plt.tight_layout(); plt.savefig('rugby_pymc.png', dpi=120, bbox_inches='tight'); plt.show()
print('max |Δattack| =', np.abs(ga-att).max().round(3), '  max |Δdefence| =', np.abs(gd-deff).max().round(3))
No description has been provided for this image 
max |Δattack| = 0.028   max |Δdefence| = 0.012

Results¶

  • PyMC (NUTS) reproduces the from-scratch sampler: home advantage ≈ 0.37 (×1.45), the same attack/defence ordering (England/Ireland strong, Scotland/Italy weak), with r̂ ≈ 1.00.
  • Both engines put every team on the y=x line in the comparison plot — the crossed-effects posterior is the same whether sampled block-by-block (Gibbs) or jointly (HMC).
  • Engine takeaway: unlike the Epil GLMM (which needed the non-centred trick to beat the funnel), here the centred sum-to-zero parameterisation samples cleanly — the team effects are well-identified by 5 matches each, so the geometry is benign.