Markov Chain Monte Carlo (MCMC): What It Is and How to Read the Diagnostics

Most Bayesian analyses cannot compute the answer directly. The posterior distribution — what you actually want, the updated belief about a parameter after seeing data — is proportional to the product of a likelihood and a prior, but turning that proportionality into an actual probability distribution requires dividing by a normalizing constant: an integral over every possible parameter value. For anything beyond the simplest models, that integral has no closed-form solution and cannot be computed exactly. Markov chain Monte Carlo (MCMC) sidesteps the problem entirely: instead of computing the posterior, it draws a large number of samples from it. Once you have enough samples, you can approximate any summary of the posterior — its mean, its spread, the probability a parameter exceeds some threshold — using ordinary sample statistics on the draws, without ever writing down the normalizing constant.

This guide explains what MCMC is doing conceptually, the algorithms researchers are most likely to encounter (Metropolis-Hastings, Gibbs sampling, and Hamiltonian Monte Carlo/NUTS), and — the part that matters most for anyone reading or reviewing a Bayesian analysis rather than writing one — how to check whether the sampler actually worked.

The core idea: sampling instead of computing

A Markov chain is a sequence of random states where each new state depends only on the current one, not on the full history before it. MCMC constructs a Markov chain specifically engineered so that its long-run stationary distribution is the posterior distribution you want to characterize. Run the chain long enough, and the sequence of parameter values it visits behaves like a set of (correlated) draws from that posterior — even though the chain never had to know the posterior’s normalizing constant to get there.

In practice this means: initialize the parameters somewhere, apply an update rule that only requires the unnormalized posterior (likelihood × prior, which you can always compute), repeat thousands of times, and treat the resulting sequence of parameter values as your sample from the posterior. Everything downstream — posterior means, credible intervals, probabilities of specific claims — is computed from that sample the same way you’d compute statistics from any dataset.

The two halves: what “Markov chain” and “Monte Carlo” each contribute

The name is a compound of two separate ideas that were developed independently and do different jobs. Reviewers and readers who find MCMC opaque usually find it opaque because the two halves have been merged before either was explained.

Monte Carlo: answer a hard integral by sampling from it

Monte Carlo on its own is the older and simpler idea: when a quantity you want is an integral or an expectation that cannot be evaluated in closed form, draw a large number of random values from the relevant distribution and compute an ordinary average over those draws instead. The average converges to the quantity you wanted as the number of draws grows. Nothing about this requires a chain, and nothing about it requires Bayesian statistics — it is a general-purpose numerical-integration trick.

The method dates to Los Alamos in the 1940s. Stanislaw Ulam traced the original idea to an intractable combinatorial calculation he attempted in 1946 — working out the probability of winning a game of solitaire — which John von Neumann took up for neutron-diffusion work, and which Nicholas Metropolis gave the name “Monte Carlo.” Metropolis and Ulam published the first paper on the method in 1949.

The limitation of plain Monte Carlo is severe and is exactly why the second half exists: it assumes you can already draw independent samples from the target distribution. For a posterior in more than a handful of dimensions, you cannot — you know the density only up to an unknown normalizing constant, which is the very thing you could not integrate.

Markov chain: get draws from a distribution you cannot sample directly

A Markov chain is a sequence of states in which the next state depends only on the current state and not on the path taken to reach it — the property usually called memorylessness. Chains with the right structure have a stationary distribution: a distribution that, once the chain is drawing from it, it keeps drawing from indefinitely.

The MCMC move is to run that logic backwards. Rather than analysing the stationary distribution of a chain you were handed, you construct a chain whose stationary distribution is the posterior you want, using an update rule that depends only on ratios of the unnormalized density. Run it long enough and the states it visits are draws from the posterior — obtained without ever computing the normalizing constant. The cost is that consecutive draws are correlated rather than independent, which is why effective sample size, not iteration count, is the number that tells you how much information the run actually contains.

The acceptance rule, concretely

In the original Metropolis form, with a symmetric proposal distribution, one iteration is: perturb the current parameter value to get a proposal; evaluate the unnormalized posterior at the proposal and at the current value. If the proposal has the higher value, accept it. If it is lower, accept it anyway with probability equal to the ratio of the two — if the posterior at the proposal is one-fifth the height of the posterior at the current value, accept with 20% probability. Otherwise the chain stays where it is and that repeated value is recorded as the next draw. Because only the ratio is ever used, the normalizing constant cancels and never has to be known.

Hastings generalized this in a 1970 Biometrika paper, defining the method for finite reversible Markov chains and giving an acceptance probability that adds a correction term for the proposal distribution, so proposals no longer have to be symmetric. That generalization is what the “Hastings” in Metropolis-Hastings refers to, and it is the reason the algorithm accommodates the asymmetric proposals real models need. Hastings also warned in the same paper that a high rejection rate signals a badly chosen proposal — the practical tuning problem samplers have wrestled with ever since.

Why the origin matters for reading the literature

The first MCMC algorithm was published by Metropolis and colleagues in the Journal of Chemical Physics in June 1953, and its target was not a posterior at all. It was the Boltzmann distribution for a system of N particles, whose normalizing constant is an integral with no closed form except in trivial cases — a statistical-physics problem with exactly the same shape as the Bayesian one. The method reached mainstream statistics only decades later, by way of Hastings (1970), Geman and Geman (1984), Tanner and Wong (1987), and above all Gelfand and Smith (1990), after which dedicated software such as BUGS appeared in the early 1990s.

The practical consequence: MCMC is a sampling algorithm for any target density known up to a constant. Bayesian posteriors are its most common target in research today, but they are not what it is for, and a paper using MCMC is not thereby a Bayesian paper.

The main algorithms

Metropolis-Hastings: propose, then accept or reject

The original and most general MCMC algorithm. At each step, propose a candidate new parameter value (often by perturbing the current value slightly), then accept or reject that candidate with a probability calculated from the ratio of the unnormalized posterior at the candidate versus the current point. Because it’s a ratio, the normalizing constant cancels out entirely — the algorithm never needs it. Accepted proposals move the chain; rejected proposals leave it in place for another iteration. Simple and general, but can be slow to explore high-dimensional or correlated parameter spaces if the proposal step size isn’t well tuned.

Gibbs sampling: update one parameter at a time

A special case that applies when you can derive the full conditional distribution of each parameter — its distribution given the current values of every other parameter. Gibbs sampling cycles through the parameters, drawing each one directly from its full conditional in turn. When those conditionals have a known, samplable form (common in many hierarchical and conjugate models), this is often more efficient than generic Metropolis-Hastings because every proposed step is accepted by construction.

Hamiltonian Monte Carlo and NUTS: using the gradient

Metropolis-Hastings and Gibbs sampling both explore the parameter space somewhat blindly, which makes them slow in the high-dimensional, correlated posteriors typical of real applied models. Hamiltonian Monte Carlo (HMC) instead treats the negative log-posterior as a physical landscape and uses its gradient to propose distant, high-probability moves — borrowing the mathematics of a particle moving through that landscape under simulated physics. This lets it explore efficiently even in many dimensions. The No-U-Turn Sampler (NUTS) is an extension that automatically tunes how far each HMC step travels, removing a parameter that otherwise has to be hand-tuned. NUTS is the default sampler in Stan and in PyMC, which is why most applied Bayesian work published today is running HMC/NUTS under the hood rather than plain Metropolis-Hastings, even when the paper just says “MCMC.”

Diagnostics: what a reviewer should actually check

A Bayesian analysis is not automatically trustworthy just because it used MCMC. The chain has to have actually converged to the target distribution and produced enough effectively independent information to summarize it — and both of those can fail silently if nobody checks. This is the section most manuscripts under-report and most reviewers don’t ask about, and it’s the part worth reading even if you never run a sampler yourself.

Warm-up / burn-in

Early in a chain’s run, its samples still reflect the (usually arbitrary) starting point rather than the target distribution. Modern software discards an initial block of iterations — called warm-up in Stan/HMC terminology or burn-in in older Metropolis/Gibbs terminology — before treating the remaining draws as posterior samples. In HMC samplers this warm-up phase also adapts the sampler’s internal tuning parameters (step size, mass matrix), so it serves double duty.

Trace plots

A trace plot shows a parameter’s sampled value across iterations, one line per chain. A well-mixing chain looks like dense, stationary noise — often described as a “fuzzy caterpillar” — with no drift, no long flat stretches, and multiple chains overlapping the same band. Trends, sticking points, or chains that never overlap are visible warning signs of a sampler that hasn’t converged.

R-hat (potential scale reduction factor)

R-hat compares the variance of a parameter’s draws within each chain to its variance across multiple independently-initialized chains. If the chains have converged to the same distribution, those two variances should agree and R-hat should sit at essentially 1.00. The original Gelman-Rubin diagnostic used a threshold of 1.1 as adequate; the widely used revised diagnostic from Vehtari, Gelman, Simpson, Carpenter, and Bürkner (Bayesian Analysis, 2021, “Rank-normalization, folding, and localization: An improved R-hat for assessing convergence of MCMC”) is stricter and also more robust to heavy-tailed or non-constant-variance posteriors that could fool the original version. The now-standard guidance is to want R-hat below 1.01 for every parameter, not just the ones a paper happens to report.

Effective sample size (ESS)

MCMC draws are autocorrelated — each sample is generated from the previous one, so consecutive draws resemble each other more than independent draws would. Effective sample size converts a chain of N correlated draws into the number of independent draws that would carry equivalent information. Four thousand correlated posterior draws might carry the statistical information of only a few hundred independent ones if the chain mixes slowly; conversely, an efficient HMC sampler can produce ESS close to the raw draw count. Vehtari et al. (2021) recommend treating a bulk-ESS (and tail-ESS, which matters for the accuracy of interval estimates) of at least roughly 400 across all chains combined as a practical minimum before trusting posterior summaries.

Autocorrelation

The underlying quantity ESS is built from: how strongly correlated a chain’s draws are with draws some number of steps earlier. High autocorrelation that decays slowly is what drives ESS down relative to the raw iteration count, and is usually a sign the sampler is taking small, inefficient steps through the posterior.

Divergent transitions (HMC/NUTS-specific)

A diagnostic unique to Hamiltonian samplers. A divergence is flagged when the simulated trajectory’s numerical integration becomes unstable, which typically happens in regions of the posterior with sharply changing curvature (a classic example is the “funnel” geometry that arises in some hierarchical models). Even a handful of divergent transitions is a signal that the reported posterior may not have been fully or accurately explored in that region — unlike R-hat and ESS, this isn’t a threshold to clear so much as a flag that should be exactly zero, or investigated and explained if it isn’t.

Multiple chains from dispersed starting points

Every diagnostic above depends on running several chains (four is a common default) initialized from different, deliberately spread-out starting values. A single chain can look perfectly well-behaved on its own trace plot while having converged to the wrong region of a multimodal posterior; only comparing independently-started chains against each other can catch that.

Priors and prior sensitivity

Every Bayesian model requires a prior distribution for each parameter, and the choice matters for what MCMC ultimately samples from. Priors are usually described on a spectrum: informative priors encode substantial existing knowledge (a previous study’s estimate, a physical constraint) and can meaningfully sharpen inference when that knowledge is well justified; weakly informative priors rule out implausible values without asserting a specific belief, and are the most common default in modern applied work; and flat (or “uninformative”) priors attempt to express no prior belief at all. Flat priors are not automatically the safe or neutral default they sound like — they can produce improper posteriors in some model structures, behave unpredictably under reparameterization, and in weakly identified models can let the prior’s implicit assumptions dominate the result anyway, just less transparently than naming a specific informative prior would. A prior sensitivity analysis — rerunning the model under a small number of alternative reasonable priors and checking whether the substantive conclusions change — is standard practice for demonstrating a result isn’t an artifact of a particular prior choice, and is increasingly expected by reviewers and journals for applied Bayesian papers.

Reporting results: credible intervals, not confidence intervals

Bayesian analyses summarize uncertainty with a credible interval: a range constructed directly from the posterior draws such that, for example, 95% of the posterior probability mass falls inside it. This has a genuinely intuitive interpretation — given the model and the prior, there is a 95% probability the true parameter value lies in that interval — which is exactly the interpretation researchers routinely (and incorrectly) apply to a frequentist confidence interval. A confidence interval instead makes a claim about the long-run behavior of the procedure across repeated sampling, not a probability statement about the specific parameter given the specific data in hand; see CASRAI’s guide to confidence intervals for the frequentist definition in full, and the guide to p-values for the parallel, equally common misinterpretation on the hypothesis-testing side. Reporting a posterior mean or median alongside a credible interval is the Bayesian equivalent of reporting a point estimate with a confidence interval or effect size in a frequentist analysis — readers need both the central estimate and a calibrated sense of its uncertainty.

Where researchers actually meet MCMC

• Hierarchical / multilevel models — MCMC is the standard fitting engine for multilevel models with partially pooled group-level parameters, including designs like cluster-randomized trials where outcomes are nested within clusters.
• Meta-analysis — Bayesian random-effects and network meta-analysis models are commonly fit with MCMC; see CASRAI’s guide to heterogeneity in meta-analysis for the frequentist I²/τ² framing these models extend.
• Item response theory and latent-variable models — Bayesian estimation is widely used for IRT models and for latent-variable structures related to confirmatory factor analysis.
• Phylogenetics — tools such as BEAST and MrBayes use MCMC to sample trees and evolutionary-rate parameters from a posterior over tree topologies, which is intractable to compute exactly for anything but the smallest datasets.
• Epidemiological and disease modelling — Bayesian MCMC is routinely used to fit transmission and other mechanistic epidemiological models to observed case data, and to propagate parameter uncertainty into downstream projections.
• Causal and regression models more broadly — Bayesian versions of the models covered in CASRAI’s regression analysis and logistic regression guides, and of causal analysis more generally, are frequently fit with MCMC when the model structure (informative priors, partial pooling, non-standard likelihoods) makes closed-form or maximum-likelihood estimation impractical.

Reproducibility: what an MCMC analysis needs to report

An MCMC analysis is not reproducible from a paper that only states “we used a Bayesian model fit with MCMC.” Because the output is a stochastic sample rather than a deterministic calculation, reproducing or auditing the result requires the paper to report, at minimum:

• The software and its version (Stan, PyMC, JAGS, and similar tools change sampler defaults and internals across versions).
• The sampler used, if not the software’s default (Metropolis-Hastings, Gibbs, HMC/NUTS).
• The number of chains, warm-up/burn-in iterations, and post-warm-up iterations retained.
• The random seed(s), or an explicit statement that results were checked for stability across seeds.
• The full prior specification for every parameter — not just the ones the paper’s headline results depend on.
• The convergence diagnostics actually obtained (R-hat, ESS, divergence count), not just a statement that “convergence was checked.”

This is the same discipline CASRAI’s guides to reproducibility infrastructure and data availability statements describe for computational work generally — a Bayesian model fit with MCMC is exactly the kind of stochastic, software-and-configuration-dependent result that needs its full specification reported to be independently checked or re-run.

Software and tools

• Stan — a probabilistic programming language with its own HMC/NUTS implementation, accessible from R, Python, and other languages.
• PyMC — a Python probabilistic programming library built on HMC/NUTS.
• brms — an R package that generates and fits Stan models from familiar regression-style formula syntax, widely used for Bayesian multilevel models without writing Stan code directly.
• JAGS — an older, still widely used engine built primarily around Gibbs sampling, accessible from R and other languages.
• NIMBLE — an R-based system for building and customizing MCMC samplers, including Gibbs, Metropolis-Hastings, and hybrid strategies.

Frequently asked questions

What problem does MCMC actually solve?

It lets researchers work with Bayesian posterior distributions that can’t be computed in closed form, by drawing samples from the posterior instead of calculating it exactly. Any summary of interest — a mean, an interval, a probability — can then be estimated from those samples.

Why do MCMC samples need “warm-up” or “burn-in” removed?

Early iterations of a chain still reflect its arbitrary starting value rather than the target posterior. Discarding an initial block of iterations before summarizing the chain avoids letting that starting-point bias contaminate the results.

What is a good R-hat value?

Under the current, more robust rank-normalized diagnostic from Vehtari et al. (2021), R-hat should be at or below roughly 1.01 for every parameter in the model, not just the headline ones. Values noticeably above that indicate the chains have not converged to the same distribution and results should not be trusted as reported.

How many MCMC chains should I run?

Four independently, randomly initialized chains is a common practical default, and is what’s needed to compute R-hat and to catch a chain that has converged to the wrong region of a multimodal posterior.

Is a 95% credible interval the same as a 95% confidence interval?

No, even though they’re often reported the same way. A 95% credible interval is a direct probability statement: given the model and prior, there is a 95% probability the parameter lies in that range. A 95% confidence interval is a statement about the procedure’s long-run behavior across repeated sampling, not a probability statement about the parameter given the data actually observed — see CASRAI’s guide to confidence intervals for the full distinction.

What does a divergent transition mean in Stan or PyMC?

It’s a warning specific to Hamiltonian/NUTS samplers indicating the numerical integration became unstable while exploring the posterior, usually because of sharply changing curvature in part of the parameter space. Divergences should be investigated (often via reparameterizing the model) rather than ignored, since they can mean part of the posterior wasn’t accurately sampled.

Are Markov chains and Monte Carlo the same thing?

No — they are two separate ideas that MCMC combines. Monte Carlo is estimating a quantity by drawing random samples and averaging over them; it assumes you can already sample from the distribution of interest. A Markov chain is a sequence of states in which each state depends only on the one before it. MCMC constructs a Markov chain whose long-run stationary distribution is the target you could not sample from directly, then applies the Monte Carlo averaging step to the states that chain visits. You can use Monte Carlo without any Markov chain (ordinary simulation from a known distribution), and Markov chains are studied extensively outside any Monte Carlo context.

Is Markov chain Monte Carlo inherently Bayesian?

No. MCMC is a general algorithm for drawing samples from any probability density that you can evaluate up to an unknown multiplicative constant. The first MCMC algorithm, published by Metropolis and colleagues in the Journal of Chemical Physics in 1953, targeted the Boltzmann distribution of an N-particle physical system, not a posterior. Bayesian posterior inference is the dominant use in research today — because posteriors are exactly densities known only up to a normalizing constant — but the method itself carries no Bayesian commitment, and MCMC appears in statistical physics, computational chemistry, and simulation-based frequentist procedures as well.

What is the Metropolis-Hastings acceptance rule?

Propose a new parameter value, then compare the unnormalized posterior at the proposal with its value at the current point. Under the original Metropolis rule (symmetric proposals), accept the proposal outright if it is higher, and otherwise accept it with probability equal to the ratio of the two heights — a proposal one-fifth as high is accepted 20% of the time. Hastings’ 1970 generalization multiplies that ratio by a correction term for the proposal distribution, which removes the requirement that proposals be symmetric. In every version the normalizing constant cancels out of the ratio, which is the whole reason the algorithm works on posteriors that cannot be normalized.

Why are MCMC draws correlated when ordinary Monte Carlo draws are not?

Because each state of the chain is generated from the previous one. A rejected proposal leaves the chain in place, and an accepted proposal is usually near where the chain already was, so adjacent draws carry overlapping information. This is not a defect — it is the price of being able to sample a distribution you could not sample directly — but it means the raw iteration count overstates how much you have learned. Effective sample size is the correct measure of information content, and it is routinely a small fraction of the number of iterations run.