Reversible-jump MCMC — inference when the model dimension is unknown

Bayesian variable selection, Part 5 — Green (1995)

Every method so far selected among models of fixed size — which coefficients, which associations. The final node asks a harder question: how many components does the model even have? When the dimension of the parameter is itself unknown, an ordinary sampler cannot move between models — the parameter vector changes length. Green's (1995) reversible-jump MCMC solves this by constructing moves that jump between dimensions while preserving detailed balance.

We reproduce Green's own flagship example: the coal-mining disaster series (191 explosions in British collieries, 1851–1962), modelled as a Poisson process whose rate is piecewise-constant with an unknown number of change-points $k$. RJMCMC infers $k$ and the change-point locations and the rates jointly. A PyMC single-switchpoint model (the special case $k=1$) cross-checks the dominant change, and a from-scratch base-R RJMCMC reproduces the posterior over $k$.

The data (rjmcmc_data.csv): year (1851–1962) and disasters (count that year).

The model, and the reversible jump

The yearly counts $y_t$ are Poisson with a step-function rate: $k$ change-points cut the span into $k+1$ segments, segment $j$ having rate $h_j$. Priors: $$ h_j \sim \text{Gamma}(a,b),\qquad k \sim \text{Poisson}(\lambda),\qquad \text{locations uniform over the between-year boundaries.} $$

Collapsing the rates. The Gamma prior is conjugate to the Poisson, so each segment's rate can be integrated out analytically, leaving a closed-form (Gamma-Poisson) marginal likelihood that depends only on the segment's disaster total $n_j$ and length $L_j$: $$ \log p(\text{segment}) = a\log b - \log\Gamma(a) + \log\Gamma(a+n_j) - (a+n_j)\log(b+L_j). $$ This is the key simplification: the between-model move then only ever changes the discrete partition, so no Jacobian is needed — the usual technical crux of RJMCMC disappears. (The rates are recovered afterwards from $h_j\mid\text{data}\sim\text{Gamma}(a+n_j,\,b+L_j)$.)

The moves. Each sweep proposes one of:

• Birth ($k\to k+1$): add a change-point at a random empty boundary — accept with $\min\!\big(1,\ e^{\Delta\mathcal L}\,\lambda/(k+1)\big)$;
• Death ($k\to k-1$): remove a random change-point — accept with $\min\!\big(1,\ e^{\Delta\mathcal L}\,k/\lambda\big)$;
• Shift: relocate a change-point (fixed dimension) — accept with $\min(1,e^{\Delta\mathcal L})$;

where $\Delta\mathcal L$ is the change in total marginal log-likelihood. The birth/death ratios are this clean because the uniform location prior and uniform proposals make the combinatorial factors cancel exactly, leaving only the likelihood ratio and the Poisson prior ratio $\lambda/(k+1)$. Hyperparameters (fixed): rate prior $a=1,\ b=1$; number-of-changes prior $\lambda=1$ (a priori "about one"). Engine: rjcp.py.

import numpy as np, pandas as pd, matplotlib.pyplot as plt
import rjcp                                                          # from-scratch reversible-jump sampler

d = pd.read_csv("rjmcmc_data.csv"); y = d["disasters"].values.astype(float); yr = d["year"].values
res = rjcp.rjmcmc(y, a=1.0, b=1.0, lam=1.0, niter=150000, burn=50000, seed=1)

kp = res["kpost"]
print("P(k >= 1) = %.3f   (a change-point is essentially certain)" % (1 - kp[0]))
print("posterior on the number of change-points k:")
for k in range(1, 5): print("  k=%d : %.3f" % (k, kp[k]))
print("\ntwo change regions: a major drop ~1890 (rate 3.1 -> ~1) and a smaller one ~1947 (-> ~0.5)")

# --- Fig: RJMCMC infers the DIMENSION (how many changes) and WHERE they are ---
fig, ax = plt.subplots(1, 2, figsize=(11, 4))

kk = np.arange(6)                                                    # (left) posterior over the NUMBER of change-points k
ax[0].bar(kk, kp[:6], color="navy", alpha=.85)
ax[0].set_xticks(kk); ax[0].grid(alpha=.25)
ax[0].set_title("Posterior over the number of change-points $k$")
ax[0].set_xlabel("k  (number of change-points)"); ax[0].set_ylabel("P(k | y)")

loc = res["loc"]                                                     # (right) posterior probability of a change-point at each year
ax[1].plot(yr, loc, color="firebrick", lw=1.5)
ax[1].fill_between(yr, loc, color="firebrick", alpha=.15)
for c in np.where(loc > 0.4 * loc.max())[0]:                         # mark the high-probability boundaries (WHERE)
    ax[1].axvline(yr[c], color="navy", ls="--", lw=1, alpha=.6)
ax[1].grid(alpha=.25)
ax[1].set_title("Change-point location posterior (where)")
ax[1].set_xlabel("year"); ax[1].set_ylabel("P(change-point at year)")

fig.suptitle("RJMCMC infers the dimension (how many) and the location (where)")
plt.tight_layout(); plt.show()

P(k >= 1) = 1.000   (a change-point is essentially certain)
posterior on the number of change-points k:
  k=1 : 0.219
  k=2 : 0.508
  k=3 : 0.204
  k=4 : 0.056

two change regions: a major drop ~1890 (rate 3.1 -> ~1) and a smaller one ~1947 (-> ~0.5)

Reading the trans-dimensional posterior

• A change-point is certain; the number is inferred, not assumed. $P(k\ge1)\approx1$ — the flat-rate model is decisively rejected — and the posterior actually favours two change-points ($k=2$: 0.51), with $k=1$ and $k=3$ also plausible (0.22, 0.20). RJMCMC delivers what a fixed-$k$ model cannot: a distribution over the number of regimes.
• Two regime shifts, not one. The location posterior (right panel) shows a major drop around 1890 — the well-known collapse in the disaster rate (from ~3.1 to ~1 per year), usually attributed to mining and safety changes — and a second, smaller decline around 1947 (down to ~0.5). The celebrated "single change-point at 1890" is really the $k{=}1$ slice of this richer picture.
• This is the whole point of reversible jump. By letting the dimension float, the sampler discovers the second regime shift that a one-switchpoint model is structurally unable to represent.

import pymc as pm                                                    # cross-check: the classic FIXED k=1 switchpoint model
T = len(y)
with pm.Model() as mod:
    tau   = pm.DiscreteUniform("tau", 0, T-1)                        # the single change-point
    early = pm.Exponential("early", 1.0); late = pm.Exponential("late", 1.0)
    rate  = pm.math.switch(np.arange(T) < tau, early, late)
    pm.Poisson("y", rate, observed=y)
    idata = pm.sample(3000, tune=2000, chains=4, random_seed=1, progressbar=False)
ty = yr[int(np.round(idata.posterior["tau"].mean().item()))]
e_ = float(idata.posterior["early"].mean()); l_ = float(idata.posterior["late"].mean())
print("PyMC single-switchpoint model (k fixed = 1):")
print("  change-point year = %d     early rate = %.2f     late rate = %.2f" % (int(ty), e_, l_))
print("-> matches the dominant change RJMCMC finds, and the classic Green (1995) result.")

# --- Fig: disaster counts with the RJMCMC posterior-mean rate profile ---
fig, ax = plt.subplots(figsize=(10, 4.5))
ax.bar(yr, y, color="firebrick", alpha=.5, width=1.0, label="disaster counts")           # data
ax.plot(yr, res["rate"], color="navy", lw=2, label="RJMCMC posterior-mean rate")         # fitted rate profile
ax.axvline(int(ty), color="navy", ls="--", lw=1.2, alpha=.8,
           label="PyMC k=1 change-point (%d)" % int(ty))
ax.grid(alpha=.25); ax.legend()
ax.set_title("Coal-mining disasters: data and RJMCMC posterior-mean rate")
ax.set_xlabel("year"); ax.set_ylabel("disasters / rate")
plt.tight_layout(); plt.show()

Multiprocess sampling (4 chains in 4 jobs)

CompoundStep

>Metropolis: [tau]

>NUTS: [early, late]

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

PyMC single-switchpoint model (k fixed = 1):
  change-point year = 1891     early rate = 3.06     late rate = 0.92
-> matches the dominant change RJMCMC finds, and the classic Green (1995) result.

Results — inference over models of different size

• The major change is unambiguous and cross-validated. Forced to use a single switchpoint, PyMC puts it at 1891 with the rate falling 3.06 → 0.92 — exactly the larger of the two shifts the reversible-jump sampler finds, and the classic Green (1995) / Carlin-Gelfand-Smith answer. The fixed-$k{=}1$ model can only find one change; RJMCMC agrees on that one and then adds the second (~1947) by letting $k$ itself be inferred.
• RJMCMC is the tool for "how many". Where SSVS and enumeration choose which terms enter a fixed-size model, reversible jump moves between models of different dimension — the number of mixture components, spline knots, or change-points. The collapsed construction here (integrate out the conjugate rates, jump over discrete partitions) sidesteps the Jacobian bookkeeping that makes RJMCMC notorious.
• Closing the arc. From spike-and-slab (Part 1) to exact enumeration (Part 2), VAR restrictions (Part 3), log-linear associations (Part 4) and now unknown dimension — the same Bayesian idea, let the data weigh the models, carried from selecting coefficients to selecting the size of the model itself.

References

• Green, P. J. (1995). Reversible jump Markov chain Monte Carlo computation and Bayesian model determination. Biometrika 82, 711–732.
• Richardson, S. & Green, P. J. (1997). On Bayesian analysis of mixtures with an unknown number of components. JRSS B 59, 731–792.
• Carlin, B. P., Gelfand, A. E. & Smith, A. F. M. (1992). Hierarchical Bayesian analysis of changepoint problems. Applied Statistics 41, 389–405.

Data: rjmcmc_data.csv — annual British coal-mining disaster counts, 1851–1962 (boot::coal). Next: rjcp_R.ipynb — the R cross-check (from-scratch base-R reversible-jump sampler).