Metropolis-Hastings: Detailed Balance and Proposal Design

A posterior density is usually available as a product, $f(\theta) = p(D \mid \theta),p(\theta)$, with the normalizing integral unavailable; that is the situation described in Conjugate Priors and the Arithmetic of Posterior Updating. Markov chain Monte Carlo answers a weaker question than an analytic posterior: it produces a sequence of states whose long-run frequency approximates the posterior, and every decision it makes compares $f$ at two points. A constant factor shared by those two values cancels in the comparison, so the unknown normalizer never has to be computed. That single observation is what makes the method applicable to the models where the closed form fails.

Why rejection and importance sampling run out

Two simpler methods also use an unnormalized density. Rejection sampling needs an envelope: a proposal distribution $g$ and a constant $M$ with $f \le M g$ everywhere, and it accepts a draw with probability $f/(Mg)$. The acceptance rate equals the target’s integral divided by $M$ times the proposal’s integral, so the rate is highest when the envelope is tight, and a tight envelope in more than a few dimensions is elusive: the envelope must cover the target’s peak while the proposal must span its tails, and the volume of the region between the two grows with dimension. Importance sampling avoids the envelope but replaces it with weights, and the weights of $n$ draws are dominated by the few draws that land where $f/g$ is largest; the effective sample size collapses as the dimension grows, and the estimator’s error is often infinite rather than merely large. MCMC replaces independent draws with dependent ones, accepting the loss of independence in exchange for a chain that adapts its location to the target as it runs.

Detailed balance produces the acceptance ratio

The construction starts from a proposition about stationary distributions. If a transition kernel $T(\theta \rightarrow \theta’)$ satisfies $\pi(\theta),T(\theta \rightarrow \theta’) = \pi(\theta’),T(\theta’ \rightarrow \theta)$ for all pairs, integrating both sides over $\theta$ shows that $\pi T = \pi$: the distribution is preserved. A kernel with that property is reversible with respect to $\pi$. Writing $T$ as a proposal step followed by an accept or reject decision, $T(\theta \rightarrow \theta’) = q(\theta’ \mid \theta),\alpha(\theta, \theta’)$, and requiring reversibility pins down the ratio of the acceptance probabilities:

$$\frac{\alpha(\theta, \theta’)}{\alpha(\theta’, \theta)} = \frac{\pi(\theta’),q(\theta \mid \theta’)}{\pi(\theta),q(\theta’ \mid \theta)}$$

Any function $\alpha$ with that ratio keeps reversibility, and the Metropolis-Hastings choice $\alpha(\theta, \theta’) = \min!\bigl(1,; r(\theta, \theta’)\bigr)$ with $r$ the right-hand side is the largest one satisfying it, so it rejects as little as the constraint allows. Since $\pi = f/Z$ and the same $Z$ appears in numerator and denominator, the target enters only as the ratio $f(\theta’)/f(\theta)$. When the proposal is symmetric, meaning $q(\theta \mid \theta’) = q(\theta’ \mid \theta)$, the proposal terms cancel entirely and the acceptance probability is $\min(1, f(\theta’)/f(\theta))$.

The implementation detail that matters is the arithmetic. Densities multiply, so their values underflow long before the interesting range of a model is reached, and a ratio of two underflowed numbers is undefined rather than small. Accepting when $\log u < \log f(\theta’) - \log f(\theta) + \log q(\theta \mid \theta’) - \log q(\theta’ \mid \theta)$ with $u$ uniform on $(0, 1]$ keeps every quantity finite and preserves the same decision.

import math, random

def log_target(x):
    # unnormalized mixture: 0.3 N(-2, 1) + 0.7 N(2, 1), log density up to a constant
    a = math.log(0.3) - 0.5 * (x + 2.0) ** 2
    b = math.log(0.7) - 0.5 * (x - 2.0) ** 2
    return max(a, b) + math.log1p(math.exp(-abs(a - b)))   # log-sum-exp of the two terms

def metropolis(steps, step_size=1.0, x=0.0):
    out = []
    for _ in range(steps):
        prop = x + random.gauss(0.0, step_size)      # symmetric proposal, so q cancels
        if math.log(random.random()) < log_target(prop) - log_target(x):
            x = prop                                 # a rejection leaves x unchanged
        out.append(x)
    return out

The proposal decides the cost per effective sample

The proposal enters the acceptance ratio even when it cancels from it: its scale sets how far each accepted step travels and how often a step is accepted, and those two quantities oppose each other. A scale so small that almost every proposal is accepted moves the chain a negligible distance per iteration, producing a sample whose autocorrelation decays slowly. A scale so large that almost every proposal is rejected leaves the chain in place, and a rejected step costs a full evaluation of the density while contributing a duplicate of the previous state. For a random-walk proposal against a high-dimensional Gaussian target, the scale that optimizes the chain’s progress has been analysed asymptotically and corresponds to an acceptance rate near 0.234; in one dimension the corresponding optimum sits noticeably higher. These are order-of-magnitude guides to a scale search, not targets to hit exactly.

Dimension changes the accounting further. A random-walk step displaces the whole parameter vector by an $O(1)$ distance while the target’s spread grows with dimension, so covering the target takes a number of steps that grows with the dimension, and if each density evaluation costs work proportional to the dimension the total cost per effective sample grows faster than linearly. Shaping the proposal like the target’s covariance removes that penalty when the covariance can be estimated, and gradient-based proposals change the scaling outright: adding a drift term proportional to $\nabla \log f$ moves the chain towards higher density, with the proposal density ratio in the acceptance rule accounting for the drift evaluated at both points. Hamiltonian dynamics with a leapfrog integrator and a Metropolis correction on the total energy accepts at a high rate for step sizes at the edge of numerical stability, at the cost of one gradient evaluation per integrator step.

Every assumption behind the guarantee

Reversibility is sufficient for stationarity and necessary for none of it, and stationarity is not convergence. The chain also has to be irreducible on the support of the target and aperiodic; periodicity occurs in practice when the proposal and the target conspire to cycle among a small set of states, and the empirical distribution of such a chain never settles. A common failure is invisible in the output: a target with well-separated modes. The chain explores the mode it starts in with correct within-mode frequencies, and a jump to another mode requires a passage through a region where $f$ is small, so the acceptance probability of such a proposal is negligible and the chain stays put. Detailed balance holds, the sample looks like a posterior, and the reported variance ignores the between-mode component entirely. Initialization interacts with the same problem: a chain started far in the tails spends its early iterations travelling towards the typical set, which is what burn-in discards, and no burn-in length derived from the target alone is available.

Adapting the proposal during the run

Because the scale that accepts and moves well is not known in advance, the temptation is to adjust it while the chain runs, raising it after a stretch of accepted steps and lowering it after rejections. Done without care this removes the property the method was built on: a chain whose transition kernel changes with its own history is not time-homogeneous, so the stationarity argument for a fixed kernel does not apply, and the sample it produces need not approach the target however long it runs. Adaptive schemes that work introduce a schedule under which the size of each adjustment shrinks as the iteration count grows; under that diminishing-adaptation condition the chain can retain the irreducibility and aperiodicity the asymptotics require, and the running averages still converge to the target. Two practical constraints follow. The adaptation must genuinely diminish rather than settle into a permanent feedback loop, and successive adjustments must be separated by long stretches of the chain, because two adjustments that both react to the same correlated sequence of acceptances push the scale in the same direction and overstate the correction. The quantity being measured is itself noisy: an acceptance rate estimated from $n$ steps has a standard error of about $\sqrt{a(1-a)/n}$ if the steps were independent, and the correlation between consecutive decisions makes the true error larger, so a rate measured over a hundred steps pins the scale to within a few percentage points at best.

The scale is not the only quantity worth choosing by measurement. A proposal that is a scalar multiple of the identity in a target with unequal variances along its axes accepts readily along the narrow directions and rarely along the wide ones; a diagonal scale vector, or a proposal shaped by an estimate of the target’s covariance, removes that mismatch. Any such estimate has to come from a preliminary run or be adapted under the diminishing schedule for the same reason as before, which leads to the usual two-stage arrangement: a short run used only to set the proposal, then a long run with the proposal frozen, from which the estimates are taken. Only the second run’s output is usable, since the first is conditioned on the very target it is exploring.

Diagnostics that separate a valid run from a plausible one

Three checks cover most of the ground. Chains started from dispersed initial values should converge to overlapping regions of the state space; disagreement between them indicates that the run has not reached the typical set, and the ratio of between-chain to within-chain variance is the usual summary of that disagreement. The autocorrelation of the recorded series gives an effective sample size, which for a series of $N$ states and autocorrelations $\rho_k$ is approximately $N / (1 + 2\sum_{k \ge 1} \rho_k)$; a scale that is too small shows up as an autocorrelation function that decays slowly. And thinning, keeping every $t$-th state, reduces storage without adding independent information, so the number of states retained should not be confused with the number of independent draws. The remaining limitation is what the method never provides: a bound on how far the sample is from the target. Where a deterministic iteration count matters more than asymptotic exactness, the optimization view in Variational Inference and the Evidence Lower Bound exchanges that guarantee for a bound of a different kind.