Markov chain Monte Carlo

Markov chain Monte Carlo constructs an ergodic Markov chain whose stationary distribution is the target law and uses its post-burn-in states as approximate dependent samples.

Reversible-jump Markov chain Monte Carlo samples a posterior on a disjoint union of model-specific parameter spaces. A dimension-matching bijection between old parameters plus auxiliary variables and new parameters plus reverse auxiliary variables supplies a reversible proposal. Its acceptance ratio includes posterior-density, model-selection, proposal-density, and Jacobian determinant factors. The equal-dimension reversible-jump acceptance probability is the ordinary density-based special case.

A pseudo-prior is a proper density assigned to a parameter that is inactive under a particular model. Integrating it out leaves that model's likelihood and marginal evidence unchanged. Its choice affects movement between model states, but not the intended marginal posterior over models.

A pseudo-prior can make model parameter vectors have equal dimensions, allowing identity-matched RJ-MCMC moves. If active and inactive coordinates share the same prior density in the two augmented targets, these factors cancel in model-switch ratios. Within-model updates must preserve the augmented target, including the inactive coordinate's pseudo-prior law.

Compare independent Poisson means with one shared mean under proper shape-rate gamma distribution priors. Adding the second mean as an inactive parameter under the shared-mean model gives a simple identity model switch. With equal model priors and identical active/inactive priors, its ratio from separate to shared means is $e^{{\gamma }_2-{\gamma }_1}({\gamma }_1/{\gamma }_2{)}^{x_2}$. Poisson-gamma conjugacy supplies exact within-model updates, and integrating the gamma kernels supplies an independent Bayes factor benchmark.

For unnormalized model posterior densities $h_j$, model-selection probabilities $s_{jk}$ and parameter-proposal densities $q_{jk}$, accept a proposal $(j,\theta )\to (k,{\theta }^{\prime })$ with the minimum of one and $h_k({\theta }^{\prime })s_{kj}({\theta }^{\prime })q_{kj}(\theta \mid {\theta }^{\prime })/[h_j(\theta )s_{jk}(\theta )q_{jk}({\theta }^{\prime }\mid \theta )]$. A deterministic matching-map formulation instead includes its explicit Jacobian determinant; it must not be counted again when already absorbed into a direct proposal density.

The hit-and-run sampler moves within a full-dimensional compact convex set by choosing a uniform direction, then a uniform point on the chord through the current point in that direction. Its invariant law is the uniform distribution on the set. Starting in the interior keeps the chain in the interior almost surely. To define absolutely continuous transitions also at boundary points, redraw directions whose chords have zero length. Without this convention, a corner can have a positive holding atom.

For a planar hit-and-run sampler, let $\ell (x,u)$ be the chord length and $a(x)$ the fraction of directions with positive chord length. With zero-length chords redrawn, the probability density function isfor almost every $y$ in the body. Interior points have $a(x)=1$. The two signed representations of a direction contribute the factor two in the polar-coordinate calculation. If the body has diameter $D$, a version of this density is bounded below by $1/(\pi D^2)$. This is a global Doeblin condition after normalization by the body's area.

Several chains from dispersed initial states, trace plots, autocorrelation, effective sample sizes, and the rank-normalized split potential scale reduction factor $\hat{R}$ assess mixing and disagreement between chains. An $\hat{R}$ close to one is useful evidence but cannot prove convergence or prove posterior propriety. Comparing chains that explore different modes is especially important; monitor each scientifically important observable and its Monte Carlo uncertainty. A long apparent plateau can also occur when an improper-target sampler is drifting toward a boundary.

Thinning keeps every $L$th state of a Markov chain. The retained chain has transition kernel $K^L$ and the same stationary distribution. It is generally still dependent, and reducing correlation between retained draws does not by itself improve precision per original transition. Keep all post-warmup draws unless storage or later processing requires thinning; assess precision using the effective sample size of a Markov chain.

A Markov kernel is geometrically ergodic when its iterates converge to stationarity in total variation at a geometric rate, with a finite state-dependent prefactor.

For an aperiodic positive Harris recurrent Markov chain with geometric ergodicity, an observable with stationary moment $\pi (|h{|}^{2+\delta })<\infty$ for some $\delta >0$ satisfies a central limit theorem. Its asymptotic variance is ${\gamma }_h(0)+2{\sum }_{k\ge 1}{\gamma }_h(k)$, with stationary covariances. Ergodicity alone does not guarantee this theorem.

Uniform geometric ergodicity strengthens geometric ergodicity by bounding the prefactor uniformly over the starting state. If $P(x,\cdot )\ge \epsilon \mu (\cdot )$ globally and $\mu$ is invariant, thenWrite $P=\epsilon \Pi +(1-\epsilon )R$, where $\Pi$ draws independently from $\mu$, to obtain the bound.

Applying a Markov kernel to two initial laws cannot increase their total variation distance, even for the distribution of an entire subsequent finite trajectory. If a bounded summand has absolute value at most $B$, the variance of a length-$N$ path sum differs between these laws by at most $6B^2{N}^2$ times their total variation distance. Thus a geometrically long enough burn-in makes this difference negligible relative to $N$.

For a reversible Markov chain with invariant law $\mu$, its Markov kernel acts as a self-adjoint operator $P$ on $L^2(\mu )$. An absolute spectral gap means that $\Vert P\Vert$ on the mean-zero functions is at most $r<1$. Geometric ergodicity implies this gap for a reversible chain. Indeed, for bounded mean-zero $g$ supported where the geometric prefactor is bounded, the total-variation estimate gives $⟨g,P^{2n}g⟩\le C_g{r}^{2n}$. The spectral theorem for normal operators puts its spectral measure in $[-r,r]$. Such test functions are dense in the mean-zero subspace.

If $\Vert P{\Vert }_{L_0^2(\mu )}\le r<1$, a stationary reversible Markov chain satisfiesConsequently the variance of a sum of $N$ consecutive values is at most $N(1+r)/(1-r)$ times the stationary variance of one value. This bounds the Markov chain Monte Carlo asymptotic variance uniformly over square-integrable functions.

For a stationary ergodic Markov chain with invariant distribution $\pi$ and an integrable function $h$, the sample mean $n^{-1}{\sum }_{t=1}^{n}h(X_t)$ converges almost surely to $\int hd\pi$. Geometric ergodicity supplies standard sufficient conditions for this conclusion and stronger limit theorems.

A minorization condition gives a common component of transition laws: $K^m(x,\cdot )\ge \alpha \nu (\cdot )$ for every $x$ in a specified set, some $m\ge 1$, $\alpha >0$, and a probability measure $\nu$.

A Doeblin condition is a minorization of the whole state space. It implies uniform geometric convergence in total variation distance.

A set $A$ is $\alpha$-small when some iterate satisfies $K^m(x,\cdot )\ge \alpha \nu (\cdot )$ for every $x\in A$ and some probability measure $\nu$.

A geometric drift $KV\le \lambda V+b{1}_A$ with $\lambda <1$, together with a small-set minorisation on $A$ and irreducibility and aperiodicity, implies geometric ergodicity.

A geometric drift condition requires a measurable $V\ge 1$, a small set $C$, constants $\lambda <1$ and $b<\infty$, andTogether with irreducibility and aperiodicity, it implies geometric ergodicity. A global Doeblin condition allows $C$ to be the entire state space and $V=1$.

Suppose a Markov kernel satisfies $PV(x)\le M<\infty$ for a measurable $V\ge 1$. For $0<\rho <1$ and $R>M/\rho$, the sublevel set $C=\{V\le R\}$ gives $PV\le \rho V+M{1}_C$. If $C$ is a small set, this is a geometric drift condition. The hypotheses of an irreducible Markov chain and an aperiodic Markov chain are also required for the usual geometric ergodicity theorem.

For a stationary chain, the asymptotic variance of the sample mean is ${\mathrm{Var}}_{\pi }\psi +2{\sum }_{k\ge 1}{\mathrm{Cov}}_{\pi }(\psi (X_0),\psi (X_k))$ when the series converges.

For a reversible kernel $K$, the inequality ${\mathrm{Var}}_{\pi }(f)\le C{\mathcal{E}}_{K^2}(f)$ is equivalent to geometric contraction of variance under iteration of $K$.

One Markov kernel dominates another in Peskun ordering when it assigns at least as much transition probability to every measurable set not containing the current state. For reversible kernels this increases Dirichlet energy and cannot decrease the spectral gap.

Hamiltonian Monte Carlo introduces an auxiliary Gaussian momentum, approximately follows Hamiltonian dynamics by a reversible volume-preserving integrator, and applies a Metropolis correction.

Surrogate Hamiltonian Monte Carlo uses an auxiliary smooth probability density function $\nu$ to generate trajectories but applies a Metropolis–Hastings acceptance probability using the intended target $\mu$. For Gaussian momenta and an involutive Metropolis proposal generated by flipped leapfrog integration, the acceptance probability is $\min (1,\mu (x^{\prime })e^{-\Vert p^{\prime }{\Vert }^2/2}/[\mu (x)e^{-\Vert p{\Vert }^2/2}])$. Only the surrogate gradient is required during the trajectory; target density evaluations are still required at the endpoints.

The Metropolis–Hastings algorithm proposes $y$ from $q(x,\cdot )$ at state $x$ and accepts it with probabilityThis acceptance rule enforces detailed balance with the target density $\pi$.

A pseudo-marginal Metropolis–Hastings algorithm replaces an intractable likelihood function by a nonnegative unbiased likelihood estimator and includes its auxiliary randomness in the state. If $\mathbb{E}[\hat{L}(\theta ,U)]=L(\theta )$ for $U\sim m_{\theta }$, the extended target is proportional to $p(\theta )\hat{L}(\theta ,u)m_{\theta }(u)$. An independence proposal for new auxiliary randomness cancels the $m_{\theta }$ factors in the Metropolis–Hastings acceptance probability. On rejection the old estimate must be retained. Integrating the extended target gives the desired marginal posterior distribution exactly.

An unbiased likelihood estimator has expectation equal to the likelihood function at each parameter. A pseudo-marginal algorithm requires this estimator to be nonnegative. Importance sampling provides $\hat{L}=M^{-1}{\sum }_j\ell (F_j)p_{\theta }(F_j)/g_{\theta }(F_j)$ for independent random variables $F_j\sim g_{\theta }$, provided the proposal covers the target support. Relative variance affects the mixing time even though unbiasedness guarantees the correct marginal target.

For target density $\pi$ and proposal density $q$, the Metropolis–Hastings acceptance probability is $1\wedge \{\pi (y)q(y,x)/(\pi (x)q(x,y))\}$.

At the current state $x$, a Metropolis–Hastings algorithm draws a candidate state $y$ from its proposal distribution $q(\cdot \mid x)$. The Metropolis–Hastings acceptance probability corrects the proposal's bias so that the target distribution remains invariant.

An involutive Metropolis proposal uses a deterministic bijection $S$ satisfying $S^2=I$. For a volume-preserving $S$ and target probability density function $\rho$, accept $S(z)$ with probability $\min (1,\rho (S(z))/\rho (z))$. The identity $\rho (z)\alpha (z)=\min (\rho (z),\rho (S(z)))$ and a change of variables under $S$ prove detailed balance. A non-unit Jacobian determinant must be included when volume is not preserved.

For a proposal $Y=AX+Z$, with isotropic $Z\sim N(0,I)$ and an independent random orthogonal matrix $A$, invariance of the distribution of $A$ under transpose makes the proposal distribution symmetric. Indeed $|y-Ax|=|x-A^{T}y|$, and averaging the isotropic Gaussian density over an inversion-invariant law yields $q(x,y)=q(y,x)$. The Metropolis–Hastings algorithm therefore accepts with the target density ratio alone.

Let $0<\beta \le 1$, let $X\in {\mathbb{R}}^p$, and proposeThe proposal distribution is $N(\sqrt{1-{\beta }^2}X,{\beta }^2{I}_p)$. If $X\sim N(0,I_p)$ independently of $Z$, then $(X,Y)$ is a jointly multivariate normal distribution invariant under exchanging $X$ and $Y$, because both random vectors have covariance matrix $I_p$ and their cross-covariance matrices are both $\sqrt{1-{\beta }^2}{I}_p$. Consequently its density $q$ satisfieswhere $\varphi$ is the standard-normal density. Thus the proposal is reversible with respect to the standard normal distribution.

The preconditioned Crank–Nicolson algorithm, or pCN algorithm, uses a Gaussian autoregressive proposal reversible with respect to a standard normal distribution in a Metropolis–Hastings algorithm. When the target density is a likelihood times the proposal's invariant Gaussian prior density, the Gaussian factors cancel from the Metropolis–Hastings acceptance probability, leaving the likelihood ratio.

Random-walk Metropolis proposes a symmetric increment from the current state and accepts a proposal $y$ from $x$ with probability $\min \{1,\pi (y)/\pi (x)\}$.

An independence Metropolis–Hastings algorithm proposes from a fixed density $q$ independent of the current state. If ${\sup }_x\pi (x)/q(x)<\infty$, it is uniformly and hence geometrically ergodic.

The effective sample size discounts a correlated chain's raw length by its integrated autocorrelation time.

For a stationary Markov chain and a scalar observable with summable autocorrelation function $\rho (k)$, the asymptotic variance of its sample mean is $\mathrm{Var}(h){\tau }_{\mathrm{int}}/n$. Thus the effective sample size of a Markov chain is approximately $n/{\tau }_{\mathrm{int}}$. Different observables can have very different autocorrelation times.

For a joint density $f_{XY}$, a two-coordinate Gibbs update draws

For Gaussian precision $\left(\begin{array}{cc}A& C\\ C& D\end{array}\right)$, a systematic update of the first coordinate followed by the second makes the centered second coordinate an AR(1) chain with coefficient $C^2/(AD)<1$. Independent conditional Gaussian noise supplies its stable innovation. This proves convergence of the Gaussian Gibbs sampler and explains slower mixing when posterior correlation has magnitude near one.

For a strictly positive posterior on finitely many binary spins, a Random-scan Gibbs sampler uses the two local posterior weights to redraw one spin. Positive full conditional distribution probabilities allow every finite sequence of spin changes, giving an irreducible Markov chain with positive self-transition probabilities.

When an observation mean depends on the neighboring spins, changing a spin affects observations centered at its neighbors. Its own observation need not involve that spin. Selecting precisely the affected factors gives a correct local full conditional distribution without recomputing the entire posterior.

Every full-conditional coordinate update satisfies detailed balance with the target distribution. Their systematic ordered composition preserves the target but need not satisfy detailed balance, since compositions of reversible kernels need not be reversible. A fixed-probability random-scan mixture remains reversible.

A blocked Gibbs sampler draws groups of coordinates jointly from their full conditional distributions. Each update preserves the joint probability distribution; dependence between blocks can still cause a long mixing time.

On a bipartite image graph, the pixels of one color are conditionally independent given the other color. Draw all black pixels from their independent full conditionals, then draw all whites conditional on the new blacks. Each color update preserves the target law, so their composition does too.

A random-scan Gibbs sampler chooses one coordinate, commonly uniformly at random, and redraws it from its complete conditional distribution while leaving all other coordinates fixed. Each coordinate update preserves the target distribution.

A single-coordinate Gibbs sampler update preserves the remaining coordinates and independently redraws the selected coordinate from its full conditional distribution. The joint measure of old and new states factors into a marginal and two identical conditional factors, making it symmetric. A fixed-probability mixture of such kernels therefore satisfies detailed balance.

A tempered Gibbs sampler changes the coordinate-selection probabilities and compensates by drawing from modified complete conditionals. Its invariant law can differ from the desired target, so expectations under the target are recovered with explicit importance weights.

The joint target density is stationary for a full Gibbs sweep, since

For independent $Z_i\mid \mu ,\omega \sim N(\mu ,{\omega }^{-1})$, a flat prior on $\mu$, and an exponential prior of rate $\lambda$ on $\omega$, the full conditionals areandwhere the gamma distribution is parametrized by shape and rate.