Pumps — R cross-check via the Negative Binomial marginal¶

Companion to pumps_python.ipynb / pumps_pymc.ipynb. Instead of MCMC, R exploits the algebra directly: marginalising the gamma rate $\theta_i$ out of the gamma–Poisson gives a Negative Binomial, so we can fit the population $(\alpha,\beta)$ with a plain MASS::glm.nb (intercept-only, exposure as offset) and then reconstruct the shrunk per-pump rates from conjugacy — empirical Bayes, no sampler.

This is the cleanest possible demonstration of the point made in the Python notebook: Pumps is the explicit-random-effects form of the Negative Binomial. The glm.nb point estimates land in the same ballpark as the Bayesian posteriors and reproduce the identical shrinkage pattern — though with only 10 pumps they are not numerically identical (see results).

The marginal¶

With $x_i\mid\theta_i\sim\text{Poisson}(\theta_i t_i)$ and $\theta_i\sim\text{Gamma}(\alpha,\beta)$, $$x_i \sim \text{NegBin}\big(\text{mean}=t_i\,\alpha/\beta,\ \ \text{size}=\alpha\big).$$ So the NB size $=\alpha$ (the gamma shape) and $\exp(\text{intercept})=\alpha/\beta$ (the population mean rate). glm.nb with offset(log t) estimates both. The per-pump empirical-Bayes rate is then the conjugate posterior mean $\hat\theta_i=(\alpha+x_i)/(\beta+t_i)$.

In [1]:
suppressMessages(library(MASS))
d <- read.csv('pumps.csv')           # x = failures, t = operating time (1000s hrs)
x <- d$x; t <- d$t

fit <- glm.nb(x ~ 1 + offset(log(t)))
alpha <- fit$theta                       # NB size = gamma shape
rate  <- exp(coef(fit)[["(Intercept)"]]) # alpha/beta = population mean rate
beta  <- alpha / rate
cat(sprintf("NB size (= alpha)            : %.3f   (Bayes posterior mean ~0.70)\n", alpha))
cat(sprintf("exp(intercept) (= alpha/beta): %.3f   (Bayes ~0.75)\n", rate))
cat(sprintf("implied beta                 : %.3f   (Bayes ~0.93)\n", beta))

NB size (= alpha)            : 0.823   (Bayes posterior mean ~0.70)
exp(intercept) (= alpha/beta): 0.652   (Bayes ~0.75)
implied beta                 : 1.262   (Bayes ~0.93)
In [2]:
# empirical-Bayes shrunk rates from conjugacy: theta_i = (alpha + x_i)/(beta + t_i)
mle <- x / t
eb  <- (alpha + x) / (beta + t)
tab <- data.frame(pump=1:10, x=x, t=t,
                  MLE_rate=round(mle,3), EB_rate=round(eb,3))
print(tab, row.names=FALSE)

png('pumps_R.png', width=620, height=460, res=110)
sz <- 0.7 + 0.02*t
plot(mle, eb, pch=19, col='seagreen', cex=sz, xlab='MLE rate  x/t',
     ylab=expression('empirical-Bayes rate  '*hat(theta)[i]),
     main='Pumps shrinkage (point size proportional to exposure)')
abline(0, 1, lty=3); abline(h=rate, col='firebrick', lty=2)
text(mle, eb, labels=1:10, pos=4, cex=0.7)
legend('topleft', legend=c('no shrinkage (y=x)', sprintf('population rate %.2f', rate)),
       lty=c(3,2), col=c('black','firebrick'), bty='n', cex=0.85)
invisible(dev.off())
cat('wrote pumps_R.png\n')

 pump  x      t MLE_rate EB_rate
    1  5  94.30    0.053   0.061
    2  1  15.70    0.064   0.107
    3  5  62.90    0.079   0.091
    4 14 126.00    0.111   0.116
    5  3   5.24    0.573   0.588
    6 19  31.40    0.605   0.607
    7  1   1.05    0.952   0.789
    8  1   1.05    0.952   0.789
    9  4   2.10    1.905   1.435
   10 22  10.50    2.095   1.940
wrote pumps_R.png

Results¶

  • glm.nb returns size = α̂ ≈ 0.82, exp(intercept) = α/β ≈ 0.65, hence β ≈ 1.26. These are the maximum-likelihood values of the marginal NB; the from-scratch Gibbs and PyMC posterior means were ≈ 0.70 / 0.75 / 0.93.
  • Why the modest gap? glm.nb is the marginal MLE with no priors, whereas the Bayesian fits carry $\alpha\sim\text{Exp}(1)$, $\beta\sim\text{Gamma}(0.1,1)$ — and with only 10 pumps the prior is not negligible. The $\alpha$ posterior is also strongly right-skewed (95% CI [0.28, 1.34]), so its posterior mean sits below the likelihood mode that glm.nb reports. Same ballpark, same model — just MLE-vs-posterior-mean at a tiny sample size.
  • The empirical-Bayes rates $(\alpha+x_i)/(\beta+t_i)$ reproduce the identical shrinkage pattern: high-exposure pumps (4, 1) near their MLEs, low-exposure pumps (7, 9) pulled toward the population. They shrink slightly harder than the full-Bayes rates because $\hat\beta$ is a touch larger.
  • The takeaway is the connection, not the third decimal: the same gamma–Poisson conjugacy that makes the from-scratch Gibbs closed-form is what gives the NB its closed-form likelihood and the rates their closed-form posterior. From-scratch Gibbs ≈ PyMC NUTS ≈ R empirical-Bayes (NB marginal) — three routes through one model.