Markov chain Monte Carlo replaces independent posterior draws by a Markov chain whose invariant law is the posterior. Correctness is about invariance and ergodic averages, not about visual randomness of the path.
Invariant distributions
The chain must preserve the posterior distribution.
Definition. A Markov kernel \(K\) has invariant distribution \(\Pi\) if \[ \int \Pi(\mathrm{d}\theta)K(\theta,A)=\Pi(A) \] for all measurable \(A\).
Proposition. If \(K\) satisfies detailed balance, \[ \Pi(\mathrm{d}\theta)K(\theta,\mathrm{d}\vartheta) =\Pi(\mathrm{d}\vartheta)K(\vartheta,\mathrm{d}\theta), \] then \(\Pi\) is invariant.
Proof. Integrate both sides of detailed balance over \(\theta\) and over \(\vartheta\in A\). The right side becomes \(\Pi(A)\), proving invariance.
Metropolis-Hastings
Metropolis-Hastings builds a reversible kernel from a proposal and an accept-reject step.
Definition. Given proposal density \(q(\theta,\vartheta)\), accept the proposed \(\vartheta\) with probability \[ \alpha(\theta,\vartheta)=1\wedge \frac{\pi(\vartheta\mid X)q(\vartheta,\theta)} {\pi(\theta\mid X)q(\theta,\vartheta)}. \]
Proof. For \(\theta\ne\vartheta\), the transition density is \(q(\theta,\vartheta)\alpha(\theta,\vartheta)\). Multiplying by \(\pi(\theta\mid X)\) gives the minimum of \(\pi(\theta)q(\theta,\vartheta)\) and \(\pi(\vartheta)q(\vartheta,\theta)\), which is symmetric in \(\theta,\vartheta\). The holding probability handles the diagonal.
Ergodic averages
After invariance, the question becomes whether time averages converge to posterior expectations.
Theorem. For a Harris ergodic Markov chain with invariant law \(\Pi\), and integrable \(g\), \[ \frac1N\sum_{s=1}^N g(\theta_s)\to\int g\,\mathrm{d}\Pi \] almost surely.
Proof. Harris ergodicity gives convergence of the transition probabilities to \(\Pi\) and recurrence strong enough for the Markov-chain ergodic theorem. Applying that theorem to the integrable function \(g\) yields almost-sure convergence of time averages to the invariant expectation. In the stationary case this is the pointwise ergodic theorem; the general Harris case follows by regeneration or small-set approximation.
Remark. Autocorrelation changes Monte Carlo variance but not the invariant target. Poor mixing is an efficiency problem unless the chain has not reached the relevant region in practical time.
Trace plots and autocorrelation
MCMC diagnostics should be interpreted as diagnostics of computation, not diagnostics of the statistical model.
Remark. Multiple chains, rank plots, effective sample size, and autocorrelation estimates all try to detect whether the simulated averages are reliable approximations to posterior integrals.
Detailed balance and invariance
Definition. A Markov kernel \(K\) is reversible with respect to \(\pi\) if \[ \pi(\mathrm{d}x)K(x,\mathrm{d}y) =\pi(\mathrm{d}y)K(y,\mathrm{d}x) \] as measures on the product space.
Proposition. Detailed balance implies invariance: \[ \int\pi(\mathrm{d}x)K(x,A)=\pi(A). \]
Proof. Integrate the detailed-balance identity over \(x\in\mathsf X\) and \(y\in A\): \[ \int_{\mathsf X}\pi(\mathrm{d}x)K(x,A) =\int_A\pi(\mathrm{d}y)K(y,\mathsf X) =\pi(A). \]
Metropolis-Hastings kernel
Theorem. Let \(q(y\mid x)\) be a proposal and define \[ \alpha(x,y) =1\wedge \frac{\pi(y)q(x\mid y)} {\pi(x)q(y\mid x)}. \] The transition that proposes \(Y\sim q(\cdot\mid x)\), accepts with probability \(\alpha(x,Y)\), and otherwise remains at \(x\), is reversible with respect to \(\pi\).
Proof. For \(x\ne y\), the accepted transition flow is \[ \pi(x)q(y\mid x)\alpha(x,y) =\min\{\pi(x)q(y\mid x), \pi(y)q(x\mid y)\}, \] which is symmetric in \(x,y\). The rejection mass lies on the diagonal and is automatically symmetric. Hence detailed balance holds.
Ergodicity and laws of large numbers
Theorem (Harris ergodic theorem). If a Markov chain is positive Harris recurrent and aperiodic with invariant probability \(\pi\), then for every \(g\in L^1(\pi)\), \[ \frac1N\sum_{s=1}^Ng(\Theta_s) \longrightarrow\int g\,\mathrm{d}\pi \quad\text{almost surely} \] from every initial state outside a \(\pi\)-null exceptional set.
Proof outline. Regeneration theory splits a positive Harris chain into iid tours between returns to a small set. The renewal-reward theorem gives the ratio of accumulated rewards to accumulated tour lengths as \(\pi g\). Aperiodicity removes cyclic subsequences, and truncation extends the result from bounded to integrable \(g\).
Invariance alone is insufficient: a kernel that never moves preserves every distribution. Irreducibility, recurrence, and adequate mixing are the properties that turn invariance into useful integration.
Markov-chain central limit theorem
Theorem. Under geometric ergodicity and a suitable moment condition, \[ \sqrt N(\bar g_N-\pi g) \Rightarrow N(0,\sigma_{\mathrm{MCMC}}^2), \] where \[ \sigma_{\mathrm{MCMC}}^2 =\gamma_0+2\sum_{k=1}^\infty\gamma_k, \qquad \gamma_k=\operatorname{Cov}_\pi\{g(\Theta_0),g(\Theta_k)\}. \]
Proof outline. Solve the Poisson equation \[ h-Kh=g-\pi g. \] Then the partial sum decomposes into a martingale \[ \sum_{s=1}^N \{h(\Theta_s)-Kh(\Theta_{s-1})\} \] plus telescoping boundary terms. Geometric ergodicity controls the solution and makes the boundary negligible. The martingale CLT gives the limit, and expanding its variance yields the autocovariance sum.
Definition. The integrated autocorrelation time and effective sample size are \[ \tau_g =1+2\sum_{k\ge1}\rho_k, \qquad N_{\mathrm{eff},g}=\frac{N}{\tau_g}, \] where \(\rho_k=\gamma_k/\gamma_0\). They depend on \(g\): one chain has no single effective sample size for every posterior feature.
Batch-means variance estimation
Proposition. Divide \(N=ab\) draws into \(a\) consecutive batches of length \(b\), with batch means \(\bar g_j\). Under a Markov-chain CLT and standard growth conditions \(a,b\to\infty\), \[ \widehat\sigma_{\mathrm{BM}}^2 =\frac{b}{a-1} \sum_{j=1}^a(\bar g_j-\bar g_N)^2 \xrightarrow{P}\sigma_{\mathrm{MCMC}}^2. \]
The Monte Carlo standard error is \(\widehat\sigma_{\mathrm{BM}}/\sqrt N\). Reporting posterior standard deviation without MCSE says nothing about numerical precision.
Gibbs sampling and blocking
Example (correlated Gaussian target). Let \[ (\theta_1,\theta_2)\sim N\left(0, \begin{bmatrix}1&\rho\\\rho&1\end{bmatrix} \right). \] The Gibbs conditionals are \[ \theta_1\mid\theta_2\sim N(\rho\theta_2,1-\rho^2), \qquad \theta_2\mid\theta_1\sim N(\rho\theta_1,1-\rho^2). \] One full sweep induces an autoregressive coefficient \(\rho^2\) for either coordinate, so mixing becomes arbitrarily slow as \(|\rho|\to1\). Drawing both coordinates jointly would be independent across iterations.
Rank-normalized split \(\widehat R\)
With multiple chains, split each chain in half, rank-normalize pooled draws, and compare between- and within-chain variances: \[ \widehat R =\sqrt{\frac{\widehat V}{W}}. \] Values appreciably above one indicate that chains have not mixed to a common distribution. Rank normalization makes the statistic more stable for heavy-tailed targets; folded ranks help detect scale disagreement. A value near one is necessary but not sufficient: all chains can be trapped in the same region.
Failure modes that diagnostics must target
- multimodality and mode isolation;
- funnels and strong posterior curvature;
- heavy tails with slow return;
- discrete bottlenecks and label switching;
- poor parameterization or proposal scale;
- high autocorrelation in the actual decision functional.
Convergence cannot be certified from finite output alone. Reliable practice combines multiple dispersed initializations, geometry-aware algorithms, function-specific ESS and MCSE, rank plots, and problem-specific invariants or simulated-data calibration.