Pumps — conjugate hierarchical Poisson (gamma–Poisson)¶
From-scratch (almost) pure-Gibbs sampler; contrast with the Epil GLMM¶
Module: pumps_gibbs.py. The classic BUGS Vol I "Pumps" example (Gaver & O'Muircheartaigh 1987): pooling failure rates across 10 pumps. It is the conjugate hierarchical Poisson — a deliberate counterpoint to the Epil model (hpois_python.ipynb), where the random effect is non-conjugate and needs Metropolis everywhere.
Model¶
10 pumps; pump $i$ records $x_i$ failures over an exposure $t_i$ (1000s of hours): $$x_i \sim \text{Poisson}(\theta_i\, t_i),\qquad \theta_i \sim \text{Gamma}(\alpha,\beta),\qquad \alpha \sim \text{Exp}(1),\ \ \beta \sim \text{Gamma}(0.1, 1.0).$$ $\theta_i$ is pump $i$'s failure rate; $t_i$ is a fixed exposure offset. The Gamma$(\alpha,\beta)$ is the population of rates — pumps with little exposure borrow strength from it.
The algorithm — why it is (almost) all closed-form¶
A Gibbs sweep cycles through three full conditionals; two are exact thanks to gamma–Poisson conjugacy.
1. Rates $\theta_i \mid x_i,t_i,\alpha,\beta$ — exact Gamma draw. $$p(\theta_i\mid\cdot)\ \propto\ \underbrace{\theta_i^{x_i}e^{-\theta_i t_i}}_{\text{Poisson}}\ \underbrace{\theta_i^{\alpha-1}e^{-\beta\theta_i}}_{\text{Gamma prior}}\ =\ \theta_i^{\alpha+x_i-1}e^{-(\beta+t_i)\theta_i}\ \equiv\ \text{Gamma}(\alpha+x_i,\ \beta+t_i).$$ The gamma is the conjugate prior for a Poisson rate: the count $x_i$ adds to the shape, the exposure $t_i$ adds to the rate. The posterior mean $\frac{\alpha+x_i}{\beta+t_i}$ is a precision-weighted blend of the population mean $\alpha/\beta$ and the pump's own MLE $x_i/t_i$ — so low-exposure pumps shrink most toward the population.
2. Scale $\beta \mid \theta,\alpha$ — exact Gamma draw. $\beta$ is the rate of the Gamma$(\alpha,\beta)$ generating the $\theta_i$: $\prod_i \text{Gamma}(\theta_i;\alpha,\beta)\propto \beta^{n\alpha}e^{-\beta\sum\theta_i}$. With the Gamma$(a_\beta,b_\beta)$ prior, $$\beta\mid\cdot\ \sim\ \text{Gamma}\!\big(a_\beta+n\alpha,\ \ b_\beta+\textstyle\sum_i\theta_i\big).$$
3. Shape $\alpha \mid \theta,\beta$ — one Metropolis step. $\alpha$ sits inside $\beta^\alpha/\Gamma(\alpha)$, so its conditional $\propto\big(\beta^\alpha/\Gamma(\alpha)\big)^n\big(\prod_i\theta_i\big)^{\alpha}e^{-\alpha}$ is non-standard → a single random-walk Metropolis update on $\log\alpha$ (with Jacobian). That's the only non-Gibbs piece.
Why this differs from the Epil GLMM (hpois_python.ipynb)¶
| Epil (hpois) | Pumps (here) | |
|---|---|---|
| random effect | $b_j\sim N(0,\sigma^2)$ on the log-rate (additive) | $\theta_i\sim\text{Gamma}(\alpha,\beta)$ on the rate (multiplicative) |
| likelihood × RE | Poisson × log-Normal → non-conjugate | Poisson × Gamma → conjugate |
| random-effect update | RW-Metropolis per unit | exact Gamma draw |
| other params | β-block Metropolis; σ² conjugate; + location move | β exact Gamma; α one Metropolis step |
| overall | full Metropolis-within-Gibbs | (almost) pure Gibbs |
The deep reason: the gamma is conjugate to the Poisson rate, but the normal is conjugate to a Gaussian mean, not to a Poisson log-rate. Pumps puts heterogeneity directly on the rate (gamma) ⇒ everything stays closed-form; Epil puts it on the log-rate (normal — natural for regression coefficients) ⇒ the log link breaks conjugacy and forces Metropolis.
Connection to the Negative Binomial. Marginalising $\theta_i$ out of the gamma–Poisson gives $x_i\sim\text{NegBin}$ — so Pumps is the explicit-random-effects form of the gamma mixing that defines the NegBin (count_reg Section 2). The same conjugacy that makes this Gibbs closed-form is what gives the NegBin its closed-form likelihood. (Epil's log-normal random effect has no closed-form marginal — there is no standard "log-normal–Poisson" distribution — which is exactly why it needs Metropolis.)
import numpy as np, matplotlib.pyplot as plt
from pumps_gibbs import pumps_data, pumps_gibbs, summary
x, t = pumps_data()
out = pumps_gibbs(x, t, R=21000, burn=1000, seed=1)
al, be = out['alpha'], out['beta']
print(f"alpha = {al.mean():.3f} 95% CI [{np.percentile(al,2.5):.2f}, {np.percentile(al,97.5):.2f}] (BUGS ~0.70)")
print(f"beta = {be.mean():.3f} 95% CI [{np.percentile(be,2.5):.2f}, {np.percentile(be,97.5):.2f}] (BUGS ~0.93)")
print(f"population mean rate alpha/beta = {(al/be).mean():.3f}; Metropolis acceptance (alpha) = {out['accept']:.2f}\n")
th = out['theta'].mean(0); mle = x / t
print(f"{'pump':>4}{'x':>4}{'t':>8}{'MLE rate':>10}{'theta (shrunk)':>16}")
for i in range(10):
print(f"{i+1:>4}{int(x[i]):>4}{t[i]:>8.2f}{mle[i]:>10.3f}{th[i]:>16.3f}")
alpha = 0.699 95% CI [0.29, 1.32] (BUGS ~0.70) beta = 0.933 95% CI [0.19, 2.26] (BUGS ~0.93) population mean rate alpha/beta = 0.937; Metropolis acceptance (alpha) = 0.58 pump x t MLE rate theta (shrunk) 1 5 94.30 0.053 0.059 2 1 15.70 0.064 0.102 3 5 62.90 0.079 0.089 4 14 126.00 0.111 0.116 5 3 5.24 0.573 0.600 6 19 31.40 0.605 0.610 7 1 1.05 0.952 0.888 8 1 1.05 0.952 0.891 9 4 2.10 1.905 1.577 10 22 10.50 2.095 1.987
fig, ax = plt.subplots(1, 2, figsize=(12, 4.4))
pop = (al/be).mean()
sizes = 30 + 8*t # point size ∝ exposure
ax[0].scatter(mle, th, s=sizes, color='steelblue', zorder=3, edgecolor='k', lw=.4)
lo, hi = 0, mle.max()*1.05; ax[0].plot([lo,hi],[lo,hi],'k:',lw=1,label='no shrinkage (y=x)')
ax[0].axhline(pop, color='firebrick', ls='--', lw=1, label=f'population rate α/β = {pop:.2f}')
for i in range(10): ax[0].annotate(str(i+1), (mle[i], th[i]), fontsize=7, xytext=(3,3), textcoords='offset points')
ax[0].set_xlabel('MLE rate $x_i/t_i$'); ax[0].set_ylabel('posterior rate $\\hat\\theta_i$')
ax[0].set_title('Shrinkage (point size ∝ exposure $t_i$)'); ax[0].legend(fontsize=8)
# posterior population rate distribution Gamma(alpha,beta) (mixed over draws) + the theta_i
g = np.linspace(0, 2.5, 300)
from scipy.stats import gamma as gammad
dens = np.mean([gammad.pdf(g, a, scale=1/b) for a, b in zip(al[::20], be[::20])], axis=0)
ax[1].plot(g, dens, 'firebrick', lw=2, label='population Gamma(α,β)')
ax[1].plot(th, np.zeros_like(th), 'o', color='steelblue', label='posterior $\\hat\\theta_i$')
ax[1].set_xlabel('failure rate θ'); ax[1].set_yticks([]); ax[1].set_title('Estimated population of rates'); ax[1].legend(fontsize=8)
plt.tight_layout(); plt.savefig('pumps_shrinkage.png', dpi=120, bbox_inches='tight'); plt.show()
Results¶
- α ≈ 0.70, β ≈ 0.93 (matching the canonical BUGS values); the population mean failure rate — the posterior mean of α/β, which is the line drawn in the figure — is ≈ 0.94 failures per 1000 h. (The ratio of the two point estimates, 0.70/0.93 ≈ 0.75, is a different, lower summary of the same quantity.) Acceptance for the single α-Metropolis step ≈ 0.58.
- Shrinkage works as it should: high-exposure pumps (4: t=126; 1: t=94) sit almost on the y=x line — their MLEs are trusted; low-exposure pumps shrink toward the population — pump 9 (t=2.1) moves from MLE 1.91 → 1.58, pump 7 (t=1.05) from 0.95 → 0.89. The smallest, noisiest pumps borrow the most strength.
- The α posterior is wide (CI [0.29, 1.32]); with only 10 pumps the shape of the rate distribution is weakly identified — but the per-pump rates are well estimated.
Takeaways¶
- Conjugacy buys a clean sampler. Placing the random effect as a Gamma on the rate makes the rates and the scale exact Gibbs draws — only the gamma shape needs a Metropolis step. Compare Epil, where a Normal on the log-rate (needed for covariates) forced Metropolis for the coefficients and every random effect.
- It is the explicit form of the Negative Binomial. Gamma–Poisson is NB marginally — Pumps and
count_reg's NegBin are two views of the same mixing, one with the rates kept explicit (for shrinkage/per-unit inference), one integrated out (for a single dispersion parameter). - Cross-checks (PyMC, R) follow in the companion notebooks.