import numpy as np, pandas as pd, matplotlib.pyplot as plt
from scipy.special import expit
from seeds_gibbs import seeds_gibbs
from betabin_gibbs import betabin_gibbs

d = pd.read_csv('baseball.csv'); r = d['r'].to_numpy(float); n = d['n'].to_numpy(float)
truth = d['remaining_avg'].to_numpy(); raw = r / n
ob = betabin_gibbs(r, n, R=30000, burn=8000, seed=1)             # conjugate beta-binomial
pbb = ob['p'].mean(0)
og = seeds_gibbs(r, n, np.ones((len(d),1)), R=30000, burn=8000, seed=1)   # logit-normal GLMM
pln = expit(og['beta'][:,0][:,None] + og['bdraws']).mean(0)
kap = np.median(ob['kappa'])
print('beta-binomial : mu %.3f   kappa(median) %.0f   shrinkage %.0f%% toward mean   (acc %.2f)'
      % (ob['mu'].mean(), kap, 100*kap/(kap+45), ob['accept']))
print('logit-normal  : pop avg %.3f   sigma %.3f' % (expit(og['beta'][:,0]).mean(), og['sigma'].mean()))
print('two engines agree: corr %.4f, max diff %.4f' % (np.corrcoef(pbb,pln)[0,1], np.abs(pbb-pln).max()))

beta-binomial : mu 0.267   kappa(median) 268   shrinkage 86% toward mean   (acc 0.33)
logit-normal  : pop avg 0.266   sigma 0.063
two engines agree: corr 0.9972, max diff 0.0199

# James-Stein validation: predict the rest-of-season average
def rmse(x): return np.sqrt(np.mean((x - truth)**2))
gm = np.full_like(truth, raw.mean())
print('RMSE predicting rest-of-season batting average:')
print('  raw 45-AB average     : %.4f' % rmse(raw))
print('  grand mean (full pool): %.4f' % rmse(gm))
print('  beta-binomial (shrunk): %.4f  (%.0f%% better than raw)' % (rmse(pbb), 100*(1-rmse(pbb)/rmse(raw))))
print('  logit-normal  (shrunk): %.4f  (%.0f%% better than raw)' % (rmse(pln), 100*(1-rmse(pln)/rmse(raw))))

RMSE predicting rest-of-season batting average:
  raw 45-AB average     : 0.0690
  grand mean (full pool): 0.0397
  beta-binomial (shrunk): 0.0384  (44% better than raw)
  logit-normal  (shrunk): 0.0390  (43% better than raw)

fig, ax = plt.subplots(1, 2, figsize=(13, 5))
# (A) shrinkage: raw -> shrunk (both models) with the truth, per player
order = np.argsort(raw); xp = np.arange(len(d)); pop = raw.mean()
for k, i in enumerate(order):
    ax[0].plot([k, k], [raw[i], pbb[i]], color='grey', lw=.8, zorder=1)
ax[0].scatter(xp, raw[order], s=35, color='darkorange', zorder=2, label='raw (45 AB)')
ax[0].scatter(xp, pbb[order], s=30, color='steelblue', zorder=3, label='beta-binomial')
ax[0].scatter(xp, pln[order], s=12, color='seagreen', zorder=4, label='logit-normal')
ax[0].scatter(xp, truth[order], marker='x', s=45, color='black', zorder=5, label='truth (rest of season)')
ax[0].axhline(pop, color='firebrick', ls='--', lw=1, label='group mean')
ax[0].set_xticks([]); ax[0].set_xlabel('player (sorted by raw average)'); ax[0].set_ylabel('batting average')
ax[0].set_title('Shrinkage: raw → shrunk lands near the truth'); ax[0].legend(fontsize=7)
# (B) prediction error: raw vs truth (scattered) and shrunk vs truth (tight on y=x)
ax[1].plot([.15,.40],[.15,.40],'k:',lw=1)
ax[1].scatter(raw, truth, s=35, color='darkorange', label=f'raw (RMSE {rmse(raw):.3f})')
ax[1].scatter(pbb, truth, s=30, color='steelblue', label=f'beta-binomial (RMSE {rmse(pbb):.3f})')
ax[1].set_xlabel('prediction'); ax[1].set_ylabel('rest-of-season truth')
ax[1].set_title('Shrunk predictions hug the truth; raw overshoots'); ax[1].legend(fontsize=8)
plt.tight_layout(); plt.savefig('baseball.png', dpi=120, bbox_inches='tight'); plt.show()

	per-player model	conjugacy	engine
logit-normal GLMM	$\text{logit}\,p_i=\mu+b_i,\ b_i\sim N(0,\sigma^2)$	non-conjugate (Metropolis)	`seeds_gibbs` (X = ones)
beta-binomial	$p_i\sim\text{Beta}(\alpha,\beta)$	conjugate: $p_i\mid r_i\sim\text{Beta}(\alpha+r_i,\beta+n_i-r_i)$	`betabin_gibbs` (new)

field	meaning
`r`, `n`	hits in the first 45 at-bats
`remaining_avg`	rest-of-season average (the "truth" to predict)

Baseball — binomial partial pooling, two ways (Efron–Morris shrinkage)¶

Logit-normal GLMM vs the conjugate beta-binomial; the James–Stein result¶

Two parameterisations of binomial shrinkage¶

The data — Efron & Morris (1975)¶

Results & interpretation¶

Takeaways¶