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 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.
