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. 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