Reversible-Jump MCMC — Inference over Model Dimension

Python · R · Part 5 of the variable-selection arc — unknown model dimension

Overview

Part 5 of the variable-selection arc — inference when the model's size is unknown. Every method so far selected among models of fixed size: which coefficients (Part 1), which models (Part 2), which VAR restrictions (Part 3), which associations (Part 4). The final question is harder: how many components does the model even have? When the parameter's dimension 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 — the bookkeeping condition that makes any MCMC sampler valid: the chain must be as likely to make a given move as to make its reverse, so that time spent in each model is proportional to that model's posterior probability rather than to how easy it is to reach. Getting it to hold across a change of dimension is the entire difficulty Green solved. This reproduces Green's 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$ — and RJMCMC infers $k$ and the change-point locations and the rates jointly.

y_t \sim \text{Poisson}(h_{s(t)}), \quad h_j\sim\text{Gamma}(a,b), \quad k\sim\text{Poisson}(\lambda); \qquad \{\text{birth}: k\!\to\!k{+}1,\ \text{death}: k\!\to\!k{-}1,\ \text{shift}\}

Jumping dimensions — collapsed rates, no Jacobian

The yearly counts are Poisson with a step-function rate: $k$ change-points cut the span into $k+1$ segments, segment $j$ having rate $h_j\sim\text{Gamma}(a,b)$ , with $k\sim\text{Poisson}(\lambda)$ and locations uniform over the between-year boundaries. Collapsing the rates is the key simplification: the Gamma prior is conjugate to the Poisson, so each segment's rate integrates 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$ . The between-model move then only ever changes the discrete partition, so no Jacobian is needed — the usual technical crux of RJMCMC disappears (rates are recovered afterward from $h_j\mid\text{data}\sim\text{Gamma}(a+n_j, b+L_j)$ ). Each sweep proposes a birth ( $k\to k+1$ , add a change-point), a death ( $k\to k-1$ , remove one), or a shift (relocate one at fixed dimension), accepted on the change in total marginal log-likelihood times the Poisson prior ratio — the combinatorial factors cancel exactly because the location prior and proposals are uniform. Hyperparameters fixed at $a=b=1$ , $\lambda=1$ . Engine: a self-contained rjcp.py.

Results

A change-point is certain; the number is inferred, not assumed. $P(k\ge 1)\approx 1$ decisively rejects the flat-rate model, 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: a major drop around 1890 (the well-known collapse from ~3.1 to ~1 disaster/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 — by letting the dimension float, the sampler discovers the second regime shift a one-switchpoint model is structurally unable to represent. Cross-validated: a PyMC single-switchpoint model (the $k=1$ special case) puts the change at 1891 with the rate falling 3.06 → 0.92 — exactly the larger of the two shifts, and the classic Green / Carlin–Gelfand–Smith answer. Unknown dimension: where SSVS and enumeration choose which terms enter a fixed-size model, reversible jump moves between models of different dimension — change-points here, and mixture components in the companion notebook — carrying the same Bayesian idea, let the data weigh the models, from selecting coefficients to selecting the size of the model itself.

Notebooks

Python Reversible-jump MCMC from scratch — birth/death/shift moves over the number of change-points with the conjugate rates collapsed (no Jacobian), plus a PyMC single-switchpoint cross-check, on the coal-mining disaster series (Python) Open → R R cross-check — a from-scratch base-R reversible-jump sampler recovering the posterior over the number of change-points (R) Open →

Downloads

rjcp.py The from-scratch reversible-jump sampler — collapsed Gamma–Poisson segment likelihoods and birth/death/shift moves over the change-point partition (new self-contained module) rjmcmc_data.csv Annual British coal-mining disaster counts, 1851–1962 (boot::coal) — Green's (1995) flagship example

References

Green, P. J. (1995). Reversible jump Markov chain Monte Carlo computation and Bayesian model determination. Biometrika 82(4), 711–732. — reversible-jump MCMC and this coal-mining change-point example
Carlin, B. P., Gelfand, A. E. & Smith, A. F. M. (1992). Hierarchical Bayesian analysis of changepoint problems. Applied Statistics 41(2), 389–405. — the single-change-point Bayesian analysis of the coal data
Richardson, S. & Green, P. J. (1997). On Bayesian analysis of mixtures with an unknown number of components. Journal of the Royal Statistical Society B 59(4), 731–792. — reversible jump over the number of mixture components (the companion case)
Chib, S. (1998). Estimation and comparison of multiple change-point models. Journal of Econometrics 86(2), 221–241. — an alternative multiple-change-point sampler
Raftery, A. E. & Akman, V. E. (1986). Bayesian analysis of a Poisson process with a change-point. Biometrika 73(1), 85–89. — Bayesian change-point inference for a Poisson process
Jarrett, R. G. (1979). A note on the intervals between coal-mining disasters. Biometrika 66(1), 191–193. — the corrected coal-mining disaster series