""" Rugby: hierarchical Poisson with CROSSED team random effects -- from-scratch Metropolis-within-Gibbs. The PyMC "A Hierarchical model for Rugby prediction" example (Baio & Blangiardo 2010), 2014 Six Nations. Richer than Epil/Pumps: every team carries TWO random effects (attack & defence) that appear, crossed, in many matches. Compare: Pumps -- one gamma RE per unit, conjugate (pumps_gibbs.py) Epil -- one normal RE per patient, non-conjugate (hpois_gibbs.py) Rugby -- two crossed normal REs per team (here) Model ----- Each match (home team h, away team a) yields two Poisson scores: home_score ~ Poisson(theta_h), log theta_h = mu + home + att[h] - def[a] away_score ~ Poisson(theta_a), log theta_a = mu + att[a] - def[h] i.e. a team's expected score rises with its own attack and the opponent's (lack of) defence; the home side gets a fixed bonus `home`. Priors: att[t] ~ N(0, sd_att^2), def[t] ~ N(0, sd_def^2) (hierarchical), mu, home ~ N(0, prior_fixed_sd^2), sd_att^2, sd_def^2 ~ Inverse-Gamma(a0, b0). Identifiability: att and def levels are each confounded with mu (att enters every event with +, def with -). Sum-to-zero is imposed by re-centring att and def each sweep and absorbing the shift into mu (exact Gibbs translation move, as in hpois_gibbs.py). Sampler (Metropolis-within-Gibbs) --------------------------------- 1. (mu, home) -- block RW-Metropolis. 2. att[t] -- per-team RW-Metropolis, vectorised over teams (curvature-scaled step). 3. def[t] -- per-team RW-Metropolis, vectorised (enters with NEGATIVE sign). + centre att, def each sweep -> mu. 4. sd_att^2, sd_def^2 -- Inverse-Gamma conjugate. """ import numpy as np import pandas as pd def rugby_data(csv='rugby.csv', year=2014): """Return the 2-events-per-match arrays: y, scorer, opp, is_home, and the team list.""" d = pd.read_csv(csv, index_col=0) d = d[d.year == year].reset_index(drop=True) teams = sorted(set(d.home_team) | set(d.away_team)) idx = {t: i for i, t in enumerate(teams)} h = d.home_team.map(idx).to_numpy(); a = d.away_team.map(idx).to_numpy() # event rows: first the home-scoring events, then the away-scoring events y = np.concatenate([d.home_score.to_numpy(), d.away_score.to_numpy()]).astype(float) scorer = np.concatenate([h, a]) # whose attack is in play opp = np.concatenate([a, h]) # whose defence is in play is_home = np.concatenate([np.ones(len(d)), np.zeros(len(d))]) return y, scorer, opp, is_home, teams def rugby_gibbs(y, scorer, opp, is_home, T, R=20000, burn=5000, seed=0, prior_fixed_sd=10.0, a0=2.0, b0=1.0, s_fix=0.06, s_att=1.3, s_def=1.3): rng = np.random.default_rng(seed) y = np.asarray(y, float); scorer = np.asarray(scorer, int); opp = np.asarray(opp, int) is_home = np.asarray(is_home, float) sy_att = np.bincount(scorer, weights=y, minlength=T) # sum y over events where t attacks sy_def = np.bincount(opp, weights=y, minlength=T) # sum y over events where t defends mu = np.log(y.mean()); home = 0.0 att = np.zeros(T); deff = np.zeros(T); sd_att = 0.3; sd_def = 0.3 def eta(): return mu + home * is_home + att[scorer] - deff[opp] keep = R - burn MU = np.zeros(keep); HOME = np.zeros(keep); ATT = np.zeros((keep, T)); DEF = np.zeros((keep, T)) SDA = np.zeros(keep); SDD = np.zeros(keep) acc_fix = 0; acc_a = 0.0; acc_d = 0.0 for r in range(R): # 1. (mu, home) block RW-Metropolis e0 = eta(); ll = np.sum(y * e0 - np.exp(e0)) - 0.5 * (mu**2 + home**2) / prior_fixed_sd**2 mp = mu + s_fix * rng.standard_normal(); hp = home + s_fix * rng.standard_normal() ep = mp + hp * is_home + att[scorer] - deff[opp] llp = np.sum(y * ep - np.exp(ep)) - 0.5 * (mp**2 + hp**2) / prior_fixed_sd**2 if np.log(rng.random()) < llp - ll: mu, home = mp, hp; acc_fix += 1 # 2. att[t] per-team RW-Metropolis (vectorised); att enters with +sign in scorer-events e0 = eta(); m = np.exp(e0) curv = np.bincount(scorer, weights=m, minlength=T) + 1.0 / sd_att**2 prop = att + s_att * rng.standard_normal(T) / np.sqrt(curv) mp_ = np.exp(e0 + (prop - att)[scorer]) d = sy_att * (prop - att) - np.bincount(scorer, weights=(mp_ - m), minlength=T) \ - 0.5 * (prop**2 - att**2) / sd_att**2 ok = np.log(rng.random(T)) < d; att = np.where(ok, prop, att); acc_a += ok.mean() c = att.mean(); att -= c; mu += c # centre att -> mu # 3. def[t] per-team RW-Metropolis (vectorised); def enters with -sign in opp-events e0 = eta(); m = np.exp(e0) curv = np.bincount(opp, weights=m, minlength=T) + 1.0 / sd_def**2 prop = deff + s_def * rng.standard_normal(T) / np.sqrt(curv) mp_ = np.exp(e0 - (prop - deff)[opp]) d = -sy_def * (prop - deff) - np.bincount(opp, weights=(mp_ - m), minlength=T) \ - 0.5 * (prop**2 - deff**2) / sd_def**2 ok = np.log(rng.random(T)) < d; deff = np.where(ok, prop, deff); acc_d += ok.mean() c = deff.mean(); deff -= c; mu -= c # centre def -> mu (def has -sign) # 4. variances (Inverse-Gamma); sum-to-zero -> (T-1) effective df sd_att = np.sqrt(1.0 / rng.gamma(a0 + (T - 1) / 2.0, 1.0 / (b0 + 0.5 * np.sum(att**2)))) sd_def = np.sqrt(1.0 / rng.gamma(a0 + (T - 1) / 2.0, 1.0 / (b0 + 0.5 * np.sum(deff**2)))) if r >= burn: i = r - burn MU[i] = mu; HOME[i] = home; ATT[i] = att; DEF[i] = deff; SDA[i] = sd_att; SDD[i] = sd_def return dict(mu=MU, home=HOME, att=ATT, deff=DEF, sd_att=SDA, sd_def=SDD, accept_fix=acc_fix / R, accept_att=acc_a / R, accept_def=acc_d / R) def summary(draws, names=None): d = draws.reshape(draws.shape[0], -1); q = d.shape[1] names = names or [f'p{j}' for j in range(q)] return [(nm, d[:, j].mean(), d[:, j].std(), *np.percentile(d[:, j], [2.5, 97.5])) for j, nm in enumerate(names)]