Markov Chain Monte Carlo

Markov Chain Monte Carlo: Concept Architecture

Orientation and learning roadmap

Markov chain Monte Carlo makes complex probability distributions computationally accessible when direct sampling or analytic integration is impractical. This page develops the target distribution, transition kernel, invariance, Metropolis–Hastings, Gibbs sampling, Hamiltonian methods, convergence, efficiency, and Monte Carlo error. It then connects diagnostics and validation to Bayesian health-economic modeling, decision uncertainty, reproducibility, and reporting.

The target distribution determines the inferential goal

MCMC is often used to explore a Bayesian posterior distribution, but the machinery applies to any target density that can be evaluated up to a normalizing constant. In Bayesian analysis, the target combines the likelihood and prior.

[ p(\boldsymbol\theta\mid\mathbf y) =\frac{p(\mathbf y\mid\boldsymbol\theta)p(\boldsymbol\theta)} {p(\mathbf y)} \propto p(\mathbf y\mid\boldsymbol\theta)p(\boldsymbol\theta). ]

The marginal likelihood (p(\mathbf y)) is often difficult to compute, but many MCMC algorithms need only ratios of the unnormalized target. This convenience does not remove the need for a proper, scientifically defensible model.

A Markov chain moves using only its current state

A Markov chain is a sequence of random states whose next-state distribution depends on the current state through a transition kernel. Conditional on the present, earlier states do not add information about the next transition.

[ P(\boldsymbol\theta_{t+1}\in A\mid \boldsymbol\theta_t,\ldots,\boldsymbol\theta_0) =P(\boldsymbol\theta_{t+1}\in A\mid\boldsymbol\theta_t). ]

The draws are therefore dependent rather than independent Monte Carlo samples. That dependence affects precision, storage, diagnostics, and the number of iterations required.

Stationarity links the chain to the target

A distribution (\pi) is stationary for transition kernel (K) if one transition preserves it. MCMC algorithms are designed so the desired target satisfies this invariance relationship.

[ \pi(A)=\int K(\boldsymbol\theta,A),\pi(d\boldsymbol\theta). ]

Stationarity alone is not enough: a chain can preserve the target while failing to explore all of it from a given start. Irreducibility, recurrence, aperiodicity, and other regularity conditions support convergence appropriate to the state space and algorithm.

Detailed balance is sufficient but not necessary

Many algorithms enforce detailed balance, which makes probability flow between any two regions symmetric under the target. This reversibility condition implies stationarity.

[ \pi(d\theta)K(\theta,d\theta') =\pi(d\theta')K(\theta',d\theta). ]

Detailed balance is not required for every valid MCMC algorithm, and nonreversible chains can also target (\pi). The essential requirement is that the implemented transition preserves the intended target and permits adequate exploration.

Metropolis–Hastings corrects a proposal mechanism

Metropolis–Hastings proposes a candidate (\theta') from a proposal density (q(\theta'\mid\theta)), then accepts or rejects it so the target remains invariant. A rejection repeats the current state, which is a valid draw and must not be discarded as if it were missing.

[ \alpha(\theta,\theta')= \min\left{1, \frac{\pi(\theta')q(\theta\mid\theta')} {\pi(\theta)q(\theta'\mid\theta)} \right}. ]

For a symmetric proposal, the proposal terms cancel. Proposal scale controls the balance between frequent small moves and infrequent larger moves, so acceptance rate alone is not a complete efficiency measure.

A two-state example makes acceptance exact

Consider a target with (\pi(A)=0.25) and (\pi(B)=0.75), and a symmetric proposal that always proposes the other state. A move from (A) to (B) is accepted with probability one, while a move from (B) to (A) is accepted with probability (0.25/0.75=1/3).

[ K= \begin{pmatrix} 0 & 1\ 1/3 & 2/3 \end{pmatrix}. ]

The distribution ((0.25,0.75)) satisfies (\pi K=\pi). The chain is irreducible and aperiodic because both states communicate and state (B) has a positive self-transition probability.

Gibbs sampling uses full conditional distributions

Gibbs sampling updates one parameter or block at a time from its conditional distribution given the remaining parameters. Every Gibbs step is accepted because the proposal is the exact full conditional.

For parameters ((\theta_1,\theta_2)), a sweep may use:

[ \theta_1^{(t+1)}\sim p(\theta_1\mid\theta_2^{(t)},\mathbf y), \qquad \theta_2^{(t+1)}\sim p(\theta_2\mid\theta_1^{(t+1)},\mathbf y). ]

Strong posterior dependence can make coordinate updates move slowly. Blocking, reparameterization, or a different sampler may improve exploration.

Hamiltonian methods use gradients to travel efficiently

Hamiltonian Monte Carlo augments parameters with momentum and simulates approximate Hamiltonian dynamics to propose distant states with high acceptance probability. Gradient information helps the sampler avoid the random-walk behavior that can limit simple proposals in high dimensions.

The numerical trajectory uses a step size and number of integration steps, followed by a Metropolis correction for discretization error. Poor geometry, extreme curvature, boundaries, or weak parameterization can still cause divergent transitions and biased exploration.

NUTS adapts trajectory length during sampling

The No-U-Turn Sampler extends Hamiltonian Monte Carlo by stopping trajectory expansion before it doubles back on itself. During warmup it typically adapts step size and a mass matrix that scales and correlates the parameter space.

Adaptation must be separated from retained sampling so the production transition is well defined. Maximum tree-depth warnings, divergences, low energy exploration, and slow effective sampling require investigation rather than automatic iteration increases.

Warmup tunes the algorithm and leaves the initial region

Early iterations may reflect arbitrary initial values and adaptation rather than the stable sampling regime. Modern software often calls this period warmup because it can include both convergence toward typical regions and tuning of algorithm parameters.

Discarding a fixed number of iterations does not prove that the remaining chain is stationary. Warmup length, initialization, adaptation settings, and retained draws should be justified through diagnostics and sensitivity checks.

Multiple chains expose failures that one chain can hide

Independent chains started from dispersed, valid initial states provide evidence about whether the sampler reaches and mixes across the same target regions. Agreement among chains is more informative than a long single trace that remains trapped in one mode.

Initial values should be overdispersed enough to challenge exploration without violating constraints or placing chains in numerically pathological regions. Identical seeds or starts can create superficially reassuring agreement.

Autocorrelation reduces the information in retained draws

Successive MCMC draws often resemble one another, so (N) retained iterations may contain much less information than (N) independent samples. The integrated autocorrelation time summarizes this loss for a particular scalar estimand.

[ \tau=1+2\sum_{k=1}^{\infty}\rho_k, \qquad ESS\approx\frac{N}{\tau}. ]

Effective sample size is quantity-specific: means, variances, tail probabilities, and nonlinear decisions can mix at different rates. Negative autocorrelation can sometimes yield an effective size greater than the nominal count under particular estimators.

A theoretical autocorrelation example quantifies efficiency

For a stationary process with autocorrelation (\rho_k=0.8^k), the integrated autocorrelation time is ((1+0.8)/(1-0.8)=9). A chain of 1,000 retained draws therefore has a theoretical effective sample size of approximately (111.11).

If the posterior standard deviation of a quantity is (2), its Monte Carlo standard error is approximately (2/\sqrt{111.11}=0.1897). More iterations reduce this simulation error, but they do not correct bias from a chain that failed to explore the target.

Monte Carlo standard error measures simulation noise

Posterior uncertainty and Monte Carlo uncertainty answer different questions. Posterior standard deviation describes uncertainty under the model, while Monte Carlo standard error describes uncertainty in a numerical summary caused by using finitely many correlated draws.

[ MCSE(\widehat\mu)\approx \frac{SD_{posterior}}{\sqrt{ESS_{mean}}}. ]

MCSE should be small relative to the reporting precision and decision-relevant differences. A narrow posterior interval can still be estimated poorly if effective sample size is inadequate.

R-hat compares within-chain and between-chain behavior

The potential scale reduction statistic (\widehat R) compares variation within chains with variation across chains after suitable splitting and modern rank-based transformations. Values materially above one indicate that chains have not mixed to a common distribution for the monitored quantity.

An (\widehat R) near one is necessary evidence, not proof of convergence. Local failures, multimodality, tail problems, and unmonitored functions of parameters can remain even when scalar summaries look acceptable.

Trace plots reveal movement that summaries conceal

A trace plot shows sampled values against iteration for every chain. Healthy traces explore the same stable region with repeated movement across its range, while trends, long plateaus, separated bands, or abrupt scale changes signal problems.

Trace plots should be examined for important transformed parameters and decision quantities as well as raw coefficients. Dense plots can hide rare transitions, so complementary rank, autocorrelation, and energy diagnostics are useful.

Divergences identify problematic Hamiltonian trajectories

A divergent transition occurs when numerical integration cannot accurately follow the posterior geometry at the chosen settings. Divergences after warmup can indicate that some regions are being sampled inaccurately and may bias posterior summaries.

Raising the target acceptance rate can reduce step size, but persistent divergences usually require reparameterization, stronger prior information, rescaling, or model revision. Deleting divergent draws is not a valid repair.

Thinning rarely creates information

Thinning retains only every (k)-th draw and discards the rest. It can reduce storage or downstream computational burden, but it usually decreases total effective information relative to keeping all available draws.

Thinning should therefore have a specific operational reason rather than being used as a convergence remedy. Better parameterization, sampling, or computation is preferable when autocorrelation is the real problem.

Reparameterization can transform difficult geometry

Equivalent model parameterizations can produce very different sampling efficiency. Hierarchical models often benefit from choosing between centered and non-centered forms according to data informativeness and group-level variation.

Scaling predictors, separating location and scale, using appropriate constrained transforms, and reducing redundant parameters can improve geometry. Reparameterization must preserve the intended model and Jacobian adjustments where required.

Priors influence both inference and computation

Priors contribute scientific information and determine posterior behavior in weakly identified regions. Unrealistically diffuse priors can create extreme scales, funnels, or near-nonidentifiability that harm both interpretation and sampling.

Prior predictive simulation checks whether the model generates plausible data before conditioning on observations. A computationally convenient prior is not automatically substantively defensible.

Posterior predictive checks assess model adequacy

Posterior predictive simulation compares replicated data generated from the fitted model with relevant features of the observed data. It can reveal systematic misfit even when MCMC diagnostics show excellent sampling of the specified posterior.

Sampling diagnostics assess computation; predictive checks assess the model-data relationship. Neither establishes causal validity, transportability, or correct handling of selection and missingness.

MCMC estimates expectations through ergodic averages

Once regularity and convergence conditions apply, posterior expectations can be estimated using averages of retained draws. The dependence changes the uncertainty of the estimate but not the basic averaging formula.

[ E_{\pi}[h(\theta)] \approx\frac{1}{N}\sum_{t=1}^{N}h(\theta^{(t)}). ]

The function (h) may represent a parameter, predicted outcome, event probability, cost, quality-adjusted life-years, or net benefit. Diagnostics and effective sample sizes should target the functions actually used in decisions.

Health-economic calibration can use posterior sampling

Bayesian calibration combines prior information with target data to infer uncertain model parameters and their dependence. MCMC is useful when the calibrated posterior is high-dimensional, nonlinear, or available only up to a normalizing constant.

The calibration likelihood must represent target uncertainty, covariance, selection, and measurement structure. Good sampler diagnostics cannot rescue an incoherent discrepancy model or incompatible calibration targets.

Hierarchical evidence models often require MCMC

Network meta-analysis, evidence synthesis, longitudinal models, missing-data models, and multilevel cost or outcome models can generate posteriors that are difficult to integrate analytically. MCMC can propagate joint parameter uncertainty through these structures.

Exchangeability assumptions, heterogeneity distributions, prior scales, treatment coding, and correlation structures remain substantive modeling choices. Posterior precision should not be mistaken for evidence quality.

MCMC and probabilistic sensitivity analysis are distinct

MCMC draws approximate a target distribution, commonly a posterior in which successive draws are correlated. Probabilistic sensitivity analysis propagates uncertain inputs through a decision model and may use independent draws from fitted parameter distributions.

Posterior MCMC draws can supply correlated inputs to probabilistic sensitivity analysis, but the two simulation layers have different purposes and diagnostics. Treating correlated MCMC draws as independent affects Monte Carlo error estimates but may still preserve joint parameter dependence if handled correctly.

Decision summaries must carry Monte Carlo error

Cost-effectiveness probabilities, expected net benefit, credible intervals, and value-of-information quantities are functions of posterior or predictive draws. Rare tail events, small incremental differences, and nested simulations may require substantially larger effective samples than posterior means.

For willingness-to-pay threshold (\lambda), incremental net monetary benefit in draw (t) is:

[ INMB^{(t)}=\lambda\Delta E^{(t)}-\Delta C^{(t)}. ]

The posterior probability of cost-effectiveness is estimated by the proportion of draws with (INMB^{(t)}>0). Its simulation error should be checked, especially near a decision probability used for communication or policy.

Nested simulation can multiply computational noise

Some value-of-information methods or individual-level models require an inner simulation for each posterior draw. Noise from the inner model can obscure variation between parameter draws and increase computation substantially.

Common random numbers, emulation, regression-based methods, analytic integration, and balanced allocation of inner and outer draws may improve efficiency. The approximation error from any shortcut should be validated against a trusted benchmark.

Random seeds support replay but not universal identity

A recorded pseudorandom seed helps reproduce a run within the same software, algorithm, version, and execution environment. It does not guarantee bitwise-identical results after changes in libraries, parallel scheduling, hardware, or random-number generators.

Reproducibility therefore requires code, data, configuration, software versions, chain-specific seeds or streams, and initialization. Independent chains should use appropriately separated random streams rather than duplicated sequences.

Excel is unsuitable as a primary MCMC engine

A spreadsheet can illustrate a two-state chain or summarize exported draws, but it lacks the diagnostic, numerical, and reproducibility controls expected for serious Bayesian computation. Volatile random functions also make exact reruns and chain management difficult.

Production MCMC should use validated statistical software with explicit algorithms, seeds, adaptation, diagnostics, and log-density evaluation. Excel can still audit posterior summaries from a fixed exported dataset if formulas and rounding are controlled.

A disciplined MCMC workflow

Reliable MCMC begins with a coherent probability model and ends with decision-relevant summaries whose computation is demonstrably precise. More iterations are useful only after structural and sampling failures are addressed. A practical sequence is:

1. Define the target. Specify data, likelihood, priors, parameters, transforms, constraints, and estimands.
2. Check the model generatively. Use prior predictive simulation and test identification on simulated data.
3. Choose the sampler. Match Gibbs, Metropolis–Hastings, Hamiltonian, or specialized methods to the model geometry.
4. Set independent chains. Use valid dispersed starts, distinct random streams, and sufficient warmup.
5. Run diagnostic quantities. Retain algorithm warnings, acceptance information, divergences, tree depth, energy, and failures.
6. Assess convergence. Examine rank-based (\widehat R), bulk and tail ESS, trace plots, autocorrelation, and important transformations.
7. Measure simulation error. Calculate MCSE for parameters, tails, predictions, and decision outputs at reported precision.
8. Revise before extending. Reparameterize, rescale, strengthen justified priors, or change algorithms before simply adding draws.
9. Validate the model. Perform posterior predictive checks, sensitivity analyses, and comparisons with analytic or simulated benchmarks.
10. Report reproducibly. Preserve code, data, versions, seeds, settings, diagnostics, retained draws, and decision calculations.

Common errors and safeguards

MCMC output can look authoritative because it contains thousands of draws and precise summaries. The most serious failures arise when a poorly explored or misspecified target is summarized as if simulation volume guaranteed validity.

• Treating retained iterations as independent observations.
• Reporting iteration count without effective sample size or Monte Carlo standard error.
• Using a single chain and declaring convergence from a stable-looking trace.
• Treating (\widehat R) near one as proof that all relevant regions were explored.
• Deleting rejected proposals or divergent transitions from the retained sample.
• Thinning automatically to address autocorrelation.
• Increasing iterations while ignoring non-identifiability, poor scaling, or pathological geometry.
• Confusing warmup with a guaranteed removal of initialization bias.
• Checking sampler diagnostics without checking model fit or prior implications.
• Running probabilistic sensitivity analysis and MCMC as if they were the same simulation task.
• Ignoring MCSE for tail probabilities, cost-effectiveness probabilities, or value-of-information estimates.
• Claiming reproducibility from a seed without recording software, algorithm, data, and configuration.

Reproducible reporting makes posterior computation auditable

A complete record identifies the probability model, likelihood, priors, parameterization, transformations, constraints, data provenance, missing-data assumptions, software and version, algorithm, adaptation settings, chain count, initialization, warmup, retained iterations, thinning, seeds, and parallelization. It should preserve warnings, rejected or divergent transitions, tree-depth behavior, acceptance diagnostics, and reasons for reruns.

Results should report rank-based (\widehat R), bulk and tail effective sample sizes, Monte Carlo standard errors, trace or rank diagnostics, posterior predictive checks, and sensitivity analyses for important quantities. Health-economic applications should additionally show how posterior draws feed costs, outcomes, net benefit, uncertainty, and decisions, with simulation error distinguished from posterior and structural uncertainty.