Markov Chain Monte Carlo (MCMC)

Markov Chain Monte Carlo — And Why It Changed Statistics Forever

Markov Chain Monte Carlo (MCMC) is the reason Bayesian statistics went from a theoretically appealing but computationally unworkable framework to a practical engine for real-world inference. The basic problem MCMC solves is this: in Bayesian analysis, you often know the shape of the posterior distribution you need — but you can’t sample from it directly, because the normalizing constant involves an integral you can’t compute. MCMC sidesteps that requirement entirely. It builds a Markov chain whose stationary distribution is your target, then runs it long enough to collect samples that behave as if they came from that distribution.

That elegant workaround, first formalized by Nicholas Metropolis and colleagues at Los Alamos National Laboratory in 1953, has become one of the most cited ideas in computational science. Bayesian inference’s role as the modern statistical backbone is only practically realizable because MCMC exists. Without it, posterior distributions for anything beyond the most trivial models would remain analytically intractable — and the explosion of Bayesian methods across machine learning, genomics, econometrics, and epidemiology would simply not have happened.

What Is MCMC, Exactly?

Markov Chain Monte Carlo is a family of algorithms for sampling from probability distributions. The name combines two ideas. Monte Carlo refers to the use of random sampling to approximate mathematical quantities — named after the Monaco casino, a metaphor for randomness. Markov chain refers to a stochastic process where the next state depends only on the current state, not on the history leading up to it. Together, they describe a procedure: run a Markov chain that wanders through a parameter space, visiting regions with probability proportional to the target distribution. After the chain has explored sufficiently, its trajectory constitutes a set of samples from that distribution.

The critical property that makes this work is ergodicity — under mild conditions, an ergodic Markov chain eventually forgets where it started and samples proportionally from its stationary distribution, regardless of initialization. That stationary distribution is designed to be exactly the posterior you want. Probability theory foundations underpin this guarantee, and understanding them makes MCMC far less mysterious. Specifically, a chain constructed to satisfy detailed balance (also called reversibility) with respect to the target distribution will have that target as its unique stationary distribution.

Why Can’t We Just Sample Directly?

This is the question that makes MCMC’s necessity click. Direct sampling methods — like inverse CDF transformation, rejection sampling, or importance sampling — work beautifully in low dimensions. But in Bayesian inference with real models, you typically have many parameters. A hierarchical model of a clinical trial might have hundreds of participant-level effects plus population-level hyperparameters. A Bayesian neural network might have thousands. In high dimensions, direct sampling fails catastrophically: most of the probability mass concentrates in a thin “typical set” that simple sampling methods can’t efficiently find.

MCMC navigates this by using local exploration — moving one step at a time, guided by the target density — rather than trying to characterize the whole space at once. Each step is computationally cheap; the full picture emerges from the accumulated trajectory. Hypothesis testing and statistical inference traditionally relied on asymptotic approximations that break down for complex models. MCMC bypasses those approximations entirely, making exact Bayesian inference feasible in practice.

The Monte Carlo Foundation: A Brief History

Monte Carlo methods predate MCMC by a few years. They were developed during the Manhattan Project at Los Alamos in the 1940s, when physicists like Stanislaw Ulam, John von Neumann, and Nicholas Metropolis needed to simulate neutron diffusion in fissile material. The technique: use random sampling to approximate integrals too complex for deterministic computation. The name was coined by Metropolis and Ulam as a code name, inspired by the Casino de Monte-Carlo — Ulam’s uncle was a compulsive gambler there.

The leap to Markov chains came in Metropolis et al.’s 1953 paper “Equation of State Calculations by Fast Computing Machines” in the Journal of Chemical Physics. The paper introduced the algorithm for sampling from the equilibrium distribution of a thermodynamic system — which is exactly the Boltzmann distribution, a probabilistic object defined up to a normalizing constant. W.K. Hastings at the University of Toronto generalized this in 1970 to arbitrary target distributions, producing the Metropolis-Hastings algorithm that remains foundational today. Understanding probability distributions in this historical context gives you a far richer appreciation of why MCMC’s acceptance criterion works.

Markov Chains: The Mathematical Engine Behind MCMC

Before diving into specific MCMC algorithms, you need to understand what a Markov chain actually is — and what properties make it useful for sampling. This section is the mathematical bedrock. Students who skip it often find themselves confused when chains misbehave during a real analysis. Understanding the theory is not optional when using MCMC; it’s what allows you to diagnose problems and fix them. Probability distributions and random variables are prerequisites for this material.

What Is a Markov Chain?

A Markov chain is a sequence of random variables X₀, X₁, X₂, … where the conditional distribution of each variable depends only on the immediately preceding value — not on the entire history. Formally, P(Xₙ₊₁ = x | Xₙ, Xₙ₋₁, …, X₀) = P(Xₙ₊₁ = x | Xₙ). This is the Markov property, often described as “memorylessness.” The chain has no recollection of where it was two steps ago; only the current state determines where it goes next.

Markov chains can be discrete (states are countable, like weather {Sunny, Rainy, Cloudy}) or continuous (states are real-valued, like parameter values in a Bayesian model). MCMC almost always operates in continuous parameter spaces — real-valued parameters with uncountably infinite possible states. The transition kernel K(x, x’) gives the probability (density) of moving from state x to state x’. Covariance and correlation structures within a Markov chain’s trajectory determine its autocorrelation — a critical concept for MCMC efficiency.

The Stationary Distribution — Why It’s the Goal

A probability distribution π is a stationary distribution (or invariant distribution) of a Markov chain if, once the chain’s distribution equals π, it remains π at all future steps. Formally, if X ~ π and X’ is drawn by one step of the chain from X, then X’ ~ π. This is the equilibrium. In MCMC, the entire algorithm is designed so that the target posterior is the stationary distribution. Running the chain long enough lets it “converge” to this equilibrium, after which samples look like draws from the posterior.

The most important tool for constructing Markov chains with a specified stationary distribution is detailed balance: π(x) K(x, x’) = π(x’) K(x’, x). This equation says the probability of being in state x and transitioning to x’ equals the probability of the reverse journey. Any Markov chain satisfying detailed balance with respect to π has π as its stationary distribution. Both Metropolis-Hastings and Gibbs sampling are constructed precisely to satisfy this condition. Probability density functions formalize the notion of π(x) in continuous parameter spaces — the density at each point, not probability of a single point.

Key Properties for MCMC Validity

Not every Markov chain with the right stationary distribution is useful for MCMC. Three additional properties are needed. Irreducibility means the chain can reach any state from any other state in a finite number of steps — it must be able to explore the entire support of the target distribution. A chain stuck in one region of the space will never correctly represent the full posterior.

Aperiodicity means the chain doesn’t cycle periodically through states — it can return to any state at any time step, not just at multiples of some period. A periodic chain won’t converge to a stationary distribution in the standard sense.

Positive recurrence means the chain will return to any region of positive probability in finite expected time. Together, irreducibility + aperiodicity + positive recurrence guarantee that the chain is ergodic — meaning it has a unique stationary distribution and time averages converge to the expectation under that distribution. This ergodic theorem is the theoretical justification for using MCMC samples to estimate posterior means, variances, and quantiles. Total probability laws connect these chain properties to the broader probabilistic framework underlying MCMC convergence guarantees.

Mixing and Convergence: The Practical Challenge

Mixing describes how quickly a Markov chain “forgets” its starting point and begins sampling from the stationary distribution. Fast mixing means the chain converges quickly; slow mixing means it takes many iterations to explore the full distribution. Poor mixing is the main practical failure mode of MCMC — it produces correlated, unrepresentative samples that can lead to completely wrong posterior inferences.

Mixing speed depends heavily on the geometry of the target distribution and the efficiency of the proposal mechanism. A target that is highly correlated between parameters (common in hierarchical models), has heavy tails, or has multiple separated modes is difficult to mix over. This is precisely the motivation for advanced algorithms like Hamiltonian Monte Carlo — they use gradient information to make intelligent proposals that traverse the typical set efficiently, dramatically improving mixing compared to naive random-walk methods. Factor analysis and dimensionality reduction techniques often reveal the correlation structure that makes MCMC mixing challenging in multivariate posteriors.

The Metropolis-Hastings Algorithm, Gibbs Sampling, and Modern Methods

The landscape of MCMC algorithms has expanded dramatically since the 1953 Metropolis paper. Each algorithm is a different strategy for constructing a Markov chain with the desired stationary distribution. Understanding their mechanisms — not just how to run them in software — is what lets you choose correctly and diagnose failures. This section covers the classical algorithms in depth, then introduces modern gradient-based methods that have become standard in contemporary practice. Regression analysis in statistical modeling provides context for why complex posteriors arise and why simple sampling methods break down.

The Metropolis-Hastings Algorithm

The Metropolis-Hastings (M-H) algorithm is the archetype. Here’s what it does, stripped to its essentials. You’re at state θ. You propose a new state θ* by drawing from a proposal distribution q(θ*|θ) — this can be any distribution you choose, typically a Gaussian centered at the current state. You compute the acceptance ratio α = [π(θ*) q(θ|θ*)] / [π(θ) q(θ*|θ)]. You accept the proposed move with probability min(1, α), or stay at θ with probability 1 − min(1, α). Repeat.

The ratio α is comparing the “desirability” of the proposed state versus the current state, corrected for any asymmetry in the proposal. Because π appears in both numerator and denominator as a ratio, the normalizing constant cancels. This is the key: you only need π up to proportionality — only the unnormalized density, which in Bayesian inference is likelihood × prior. Bayes’ theorem directly yields this unnormalized posterior, making M-H immediately applicable to any Bayesian model where you can evaluate the likelihood and prior. The original paper by Metropolis et al. (1953) in the Journal of Chemical Physics remains the primary citation for this foundational result.

The Random Walk Metropolis: The Simplest Case

The most common special case is the random walk Metropolis, where the proposal is symmetric: q(θ*|θ) = q(θ|θ*). This causes the proposal correction terms to cancel, simplifying the acceptance ratio to just π(θ*)/π(θ) — the ratio of unnormalized target densities at the proposed and current states. Geometrically, this is intuitive: proposed moves to higher-density regions are always accepted; proposed moves to lower-density regions are accepted with probability equal to the density ratio. The chain climbs toward the posterior mode but also explores the tails probabilistically.

The critical tuning parameter for random walk Metropolis is the proposal variance. Too small, and the chain takes tiny steps — it explores slowly, samples are highly autocorrelated, and effective sample size is tiny. Too large, and most proposals land in very low-density regions and get rejected — the chain barely moves. The optimal acceptance rate (the fraction of proposals accepted) for a Gaussian random walk in high dimensions is approximately 23.4%, a theoretical result from Roberts, Gelman, and Gilks (1997). This target acceptance rate is a key diagnostic: if you’re accepting far more or far fewer moves, adjust the proposal variance. Statistical misuse in MCMC analyses often traces back to poorly tuned proposals and unchecked convergence.

# Simple Random Walk Metropolis in Python
import numpy as np

def metropolis_rw(log_target, x0, proposal_sd, n_samples):
    samples = [x0]
    x = x0
    n_accept = 0
    
    for _ in range(n_samples):
        # Propose new state
        x_proposed = x + np.random.normal(0, proposal_sd)
        
        # Acceptance ratio (log scale for numerical stability)
        log_alpha = log_target(x_proposed) - log_target(x)
        
        # Accept or reject
        if np.log(np.random.uniform()) < log_alpha:
            x = x_proposed
            n_accept += 1
        
        samples.append(x)
    
    acceptance_rate = n_accept / n_samples
    return np.array(samples), acceptance_rate

Gibbs Sampling: Coordinate-Wise Updating

Gibbs sampling takes a completely different approach. Instead of proposing a new joint state and accepting or rejecting it, Gibbs sampling updates one parameter at a time, always drawing from the exact full conditional distribution — the distribution of one parameter given all others are fixed at their current values. No acceptance-rejection step is needed; every draw is accepted automatically.

The Gibbs sampler was developed by Stuart and Donald Geman in their 1984 paper on image reconstruction using Gibbs distributions (published in the IEEE Transactions on Pattern Analysis and Machine Intelligence). Alan Gelfand and Adrian Smith at the University of Nottingham brought Gibbs sampling into mainstream Bayesian statistics in their 1990 paper in the Journal of the American Statistical Association. Gelfand and Smith (1990) is widely credited with launching the modern Bayesian computing revolution. Statistics assignment guidance on Gibbs sampling almost always cites this paper as the foundational applied reference.

How Gibbs Sampling Works Step by Step

Suppose your model has parameters θ = (θ₁, θ₂, θ₃). Starting from some initial values, one Gibbs iteration proceeds as follows: draw θ₁ from its full conditional P(θ₁ | θ₂, θ₃, data); then draw θ₂ from P(θ₂ | θ₁_new, θ₃, data); then draw θ₃ from P(θ₃ | θ₁_new, θ₂_new, data). This constitutes one complete “sweep.” The resulting sequence is a Markov chain (each state depends only on the previous) that can be shown to satisfy detailed balance with respect to the joint posterior. After burn-in, the samples represent draws from the joint posterior of all parameters simultaneously.

The power of Gibbs sampling in practice comes from conjugate priors — prior distributions that, combined with specific likelihood functions, produce full conditionals from the same distributional family. The normal-normal, beta-binomial, Dirichlet-multinomial, and gamma-Poisson conjugate pairs all yield analytically tractable full conditionals that can be sampled directly. JAGS (Just Another Gibbs Sampler), developed by Martyn Plummer at the International Agency for Research on Cancer, automates this process — it parses a model specification and automatically constructs a Gibbs sampler where possible, falling back to Metropolis steps where full conditionals lack closed forms. The beta distribution and gamma distribution are among the most important conjugate families in this framework.

Hamiltonian Monte Carlo (HMC): The Modern Standard

Hamiltonian Monte Carlo (HMC) is a fundamentally different approach that exploits gradient information to make efficient, large-step proposals. The key problem with random walk algorithms in high dimensions is that they explore by diffusion — random bumbling that scales very poorly as the number of parameters increases. The typical set becomes increasingly thin and hard to find. HMC addresses this by borrowing mechanics from physics.

The idea: treat the parameter vector θ as a position in physical space, introduce an auxiliary momentum variable p of the same dimension, and define a total energy function (Hamiltonian) H(θ, p) = -log π(θ) + p²/2. Then simulate the Hamiltonian dynamics — the path a hypothetical particle would follow under this energy — for some number of steps using the leapfrog integrator. The resulting proposed state is far from the starting point but has essentially the same total energy, making acceptance highly likely. The gradient of the log-posterior guides the trajectory away from low-density regions. Causal inference models in economics and epidemiology have benefited enormously from HMC’s ability to handle their complex posterior geometries.

The No-U-Turn Sampler (NUTS) — Stan’s Algorithm

HMC’s main tuning challenge is choosing the trajectory length — how many leapfrog steps to take before making the proposal. Too few and it behaves like a random walk; too many and the trajectory doubles back (“U-turns”), wasting computation. Matthew Hoffman and Andrew Gelman at Columbia University solved this with the No-U-Turn Sampler (NUTS), published in the Journal of Machine Learning Research in 2014. NUTS automatically identifies the point where the trajectory would start turning back, terminating there. This eliminates the trajectory length as a manual tuning parameter, making HMC fully automatic. Hoffman and Gelman (2014) is the definitive reference for this algorithm, now standard in Stan and PyMC.

The performance difference is dramatic. In models with 100+ parameters, well-tuned HMC/NUTS can produce effective sample sizes 10–100× larger than random walk Metropolis for the same number of iterations. Multiple regression models with many correlated predictors — which create challenging posterior geometries — are a prime case where HMC massively outperforms classical MCMC methods.

MCMC and Bayesian Inference: How They Work Together

Markov Chain Monte Carlo and Bayesian inference are inseparable in practice. Bayesian methods offer a principled framework for incorporating prior knowledge and quantifying uncertainty — but they produce posterior distributions that are often analytically intractable. MCMC is the computational tool that makes Bayesian inference feasible for real models. Understanding this relationship at a deep level is essential for any serious statistical analysis. Bayesian inference as the backbone of modern statistics is realized specifically through MCMC’s computational power.

Bayes’ Theorem and the Intractability Problem

Bayes’ theorem states that the posterior distribution of parameters θ given data y is: P(θ|y) = P(y|θ) P(θ) / P(y). The numerator — likelihood times prior — is straightforward to evaluate for most models. The denominator, P(y) = ∫ P(y|θ) P(θ) dθ, is the marginal likelihood. This integral is over all possible parameter values, which in high dimensions is computationally intractable except for models with conjugate priors or special structure.

MCMC bypasses P(y) entirely. Since it’s a constant with respect to θ, the acceptance ratio in Metropolis-Hastings becomes [P(y|θ*) P(θ*)] / [P(y|θ) P(θ)] — P(y) cancels. The chain samples from P(θ|y) without ever computing the denominator. This is profound: Bayesian inference for arbitrarily complex models becomes computable, limited only by the speed of evaluating the likelihood at each proposed parameter value. Applications of Bayes’ theorem in modern science are almost exclusively implemented via MCMC or related methods like variational inference.

Prior Distributions and Their Effect on MCMC

In Bayesian modeling, the prior distribution P(θ) encodes beliefs about parameters before seeing data. Priors influence MCMC behavior in important ways beyond just their effect on the posterior. Weakly informative priors — like the normal and half-normal priors recommended by Andrew Gelman and the Stan development team at Columbia University — constrain parameters to reasonable ranges without strongly determining the result. They help MCMC mixing by preventing the chain from exploring implausible parameter regions where the likelihood is near zero.

Improper priors (unnormalized, like a flat uniform prior over all reals) can create proper posteriors when combined with informative data, but they can also produce improper posteriors that MCMC chains will fail to represent correctly — often manifesting as chains that drift without bound. Conjugate priors — where the prior and posterior are from the same family — are especially useful for Gibbs sampling because they guarantee analytically tractable full conditionals. The Stan documentation from the Stan Development Team provides comprehensive guidance on prior selection for efficient MCMC sampling. Beta distributions as conjugate priors for proportions and gamma distributions as conjugate priors for rates are among the most important examples in applied Bayesian modeling.

Hierarchical Models: Where MCMC Shines

Hierarchical Bayesian models (also called multilevel models) are where MCMC demonstrates its greatest advantage over frequentist alternatives. In a hierarchical model, individual-level parameters (e.g., student performance in each school) are drawn from a group-level distribution (school effects), which are in turn drawn from a population-level distribution (educational system). This structure naturally represents partial pooling — sharing statistical strength across groups.

The posterior in hierarchical models has complex structure: the individual and group parameters are often highly correlated, creating the geometry that challenges random walk samplers. Andrew Gelman and colleagues at Columbia University developed a reparametrization called non-centered parameterization that transforms the model so these correlations are reduced, dramatically improving HMC/NUTS mixing. This is implemented automatically in Stan and discussed extensively in Betancourt and Girolami (2015). Multivariate analysis in educational and social science research routinely deploys hierarchical MCMC models for precisely this reason.

Posterior Predictive Checking: Using MCMC Samples for Model Evaluation

Once you have MCMC samples from the posterior, you can use them for posterior predictive checking — simulating new datasets from the fitted model and comparing them to the observed data. If the model is well-specified, simulated data should look statistically similar to observed data. Systematic discrepancies reveal model misspecification. This is one of the most powerful uses of MCMC samples beyond point estimation: because you have the full posterior rather than just a point estimate, you can propagate uncertainty through any downstream quantity. Residual analysis in statistical modeling parallels posterior predictive checking in purpose — both assess how well the model fits the data.

MCMC Convergence Diagnostics: How to Know Your Chain Has Worked

Running MCMC is easy. Knowing whether it worked is harder. This is where many students and even experienced practitioners make critical errors — they run a chain, get numbers out, and use those numbers for inference without confirming the chain has actually converged to the target distribution. The consequences of using non-converged chains are severe: posterior summaries can be completely wrong, systematically biased by the chain’s starting point or trapped in a non-representative region of the space. Type I and Type II error considerations extend to MCMC — invalid inference from non-converged chains is a form of systematic error that conventional frequentist correction doesn’t address.

Why Convergence Assessment Is Not Optional

Unlike optimization algorithms, MCMC algorithms don’t “converge” in the sense of reaching a fixed point. The chain never stops moving — it keeps wandering through the posterior. “Convergence” in MCMC means the chain’s distribution has stabilized to match the target distribution. There’s no moment where an alarm goes off to say “done.” You must actively check, using multiple diagnostics in combination. No single diagnostic is sufficient, and all diagnostics can miss problems in edge cases. Critical thinking in statistical assignments means applying skepticism to your own results — especially when they’re suspiciously clean.

The Gelman-Rubin R-hat Statistic

The R-hat statistic (also written R̂ or potential scale reduction factor, PSRF) was developed by Andrew Gelman and Donald Rubin at Harvard University, published in Statistical Science in 1992. It compares the variance of samples within each chain to the variance between chains when multiple chains are run from different starting points. The logic: if all chains have converged to the same distribution, within-chain variance should approximately equal between-chain variance, and R-hat should be close to 1.0.

The traditional threshold is R-hat < 1.1, but contemporary best practice (as recommended by Aki Vehtari and colleagues in their 2021 update published in Bayesian Analysis) uses R-hat < 1.05 for greater stringency, alongside bulk-ESS and tail-ESS measures. Values significantly above 1.1 indicate the chains have not mixed — they’re exploring different regions of the posterior and haven’t converged. Descriptive and inferential statistics principles apply directly to the sample summaries used in R-hat calculation.

Trace Plots

Trace plots are the most intuitive convergence diagnostic — they plot the value of each parameter against iteration number for each chain. A well-converged, well-mixed chain looks like a “fuzzy caterpillar”: dense, roughly stationary, with all chains overlapping and fluctuating around the same central value. Problematic traces reveal specific failure modes: a chain that drifts steadily upward or downward (trending — not yet converged); chains that separate and explore different regions without overlapping (multi-modal posterior or poor mixing); a chain with long “runs” where it stays near the same value for many steps (high autocorrelation, poor mixing, possibly stuck near a posterior boundary).

Stan, PyMC, and ArviZ (the Bayesian visualization library for Python) all produce trace plots automatically. In R, the bayesplot package and the coda package provide diagnostic plotting tools. Visual inspection of trace plots remains irreplaceable — automated diagnostics catch many problems, but an experienced eye catches nuances that numbers miss. Data visualization techniques for statistical analysis apply directly to interpreting MCMC trace plots and posterior distribution histograms.

Effective Sample Size (ESS)

Because MCMC samples are autocorrelated (successive samples are not independent), having 10,000 MCMC draws does not mean you have 10,000 independent pieces of information. The effective sample size (ESS) quantifies how many independent samples the correlated chain is equivalent to. If ESS = 500 from 10,000 draws, the autocorrelation has consumed 95% of the information — you effectively have 500 independent samples, not 10,000.

Low ESS manifests in imprecise posterior estimates and wide Monte Carlo error bands. The rule of thumb: you need ESS ≥ 400 for reliable bulk posterior inference, and ESS ≥ 400 specifically in the tails for reliable tail probability estimation (credible interval endpoints). Modern ESS calculations distinguish bulk-ESS (central tendency) from tail-ESS (quantiles), as autocorrelation often affects the tails differently from the bulk. Cross-validation and bootstrapping methods for assessing sampling uncertainty are conceptually related to ESS in that both address how much information a sample truly contains. Stan reports both ESS metrics automatically, flagging values below 400 as warnings.

Autocorrelation Plots and Thinning

Autocorrelation plots show the correlation between samples separated by increasing lags. At lag 1, samples are most correlated; the autocorrelation should decay toward zero as lag increases. Fast decay (by lag 10–20) indicates good mixing; slow decay (autocorrelation still significant at lag 100+) indicates poor mixing and low ESS. Thinning — keeping only every kth sample — can reduce storage and apparent autocorrelation, but it doesn’t increase information: the ESS of the thinned chain is the same as the unthinned chain. The modern consensus, championed by Gelman and colleagues, is that thinning is usually wasteful — keeping all samples and computing ESS is better than discarding samples to reduce apparent autocorrelation.

Diagnosing Specific Failure Modes

Different convergence failures have different signatures. Divergent transitions — unique to HMC/NUTS — occur when the numerical integrator encounters regions of very high curvature in the log-posterior, causing the trajectory to diverge from the true Hamiltonian dynamics. Even a few divergent transitions indicate a posterior geometry problem (often high correlation or heavy tails) that makes inference unreliable. Stan reports divergences explicitly and recommends reparametrization or adjusting the target acceptance probability. Model selection criteria like AIC and BIC become relevant after convergence — if your model doesn’t converge well, selection among well-converged alternatives is the correct next step.

Multi-modal posteriors are perhaps the most dangerous failure mode: a chain may appear to have converged perfectly within one mode while completely missing another mode with substantial posterior probability. Multiple chains from different starting points and careful examination of the full joint posterior can help detect this, but it’s never fully guaranteed. Tempering methods (parallel tempering, simulated tempering) can help chains jump between modes, but they add significant computational complexity. Survival analysis and complex model fitting frequently encounters multimodal posteriors when mixture components or competing parametrizations are possible.

MCMC Software: Stan, PyMC, JAGS, and the Modern Ecosystem

Understanding MCMC algorithms theoretically is one thing; implementing them effectively requires knowledge of the software ecosystem. Modern probabilistic programming languages abstract away the low-level MCMC machinery, letting you focus on model specification while the software handles sampling. But these tools make choices that you need to understand to use them correctly. Machine learning and statistical modeling in Python and R increasingly rely on this probabilistic programming ecosystem, making fluency with these tools an important professional skill.

Stan: The Current Gold Standard

Stan is an open-source probabilistic programming language and Bayesian inference platform developed by the Stan Development Team, originally at Columbia University by Andrew Gelman, Bob Carpenter, Matt Hoffman, and colleagues. Stan uses Hamiltonian Monte Carlo with the NUTS sampler as its default algorithm — the same algorithm described in Section 3. It performs automatic differentiation to compute gradients of the log-posterior, which HMC requires. Stan interfaces are available for R (RStan), Python (PyStan), Julia, Stata, and MATLAB. The CmdStan interface provides a command-line interface for maximum performance.

What makes Stan uniquely powerful is its combination of efficient sampling (HMC/NUTS), comprehensive diagnostics (R-hat, ESS, divergences, energy diagnostics), and a rich modeling language that handles continuous parameters, transformed parameters, generated quantities, and hierarchical structures naturally. The Stan User’s Guide is exceptionally well-written and is both a software reference and a practical Bayesian modeling textbook. Statistics assignment help for Stan models requires understanding both the modeling language and the MCMC backend.

// Stan model: Normal linear regression
data {
  int<lower=0> N;           // number of observations
  vector[N] x;              // predictor
  vector[N] y;              // outcome
}

parameters {
  real alpha;               // intercept
  real beta;                // slope
  real<lower=0> sigma;      // error SD
}

model {
  // Weakly informative priors
  alpha ~ normal(0, 10);
  beta ~ normal(0, 10);
  sigma ~ half_normal(0, 1);
  
  // Likelihood
  y ~ normal(alpha + beta * x, sigma);
}

PyMC: Bayesian Inference in Python

PyMC (formerly PyMC3) is a Python library for Bayesian statistical modeling, originally developed by Chris Fonnesbeck at Vanderbilt University. It provides an intuitive Python-native API for defining probabilistic models and supports multiple inference backends: NUTS (via Aesara or JAX), ADVI (automatic differentiation variational inference), and SMC (Sequential Monte Carlo). PyMC is the most popular choice in the Python data science ecosystem and integrates seamlessly with NumPy, pandas, and ArviZ for diagnostics and visualization.

The current version (PyMC v5+) uses the JAX backend for computation, enabling GPU acceleration — dramatically faster MCMC for large models. ArviZ, developed as the companion visualization library, produces publication-quality posterior plots, trace plots, forest plots, and diagnostic summaries. Together, PyMC + ArviZ constitute a complete Bayesian workflow environment in Python. Data science assignments using Bayesian methods almost always use PyMC as the implementation framework.

import pymc as pm
import numpy as np

# PyMC normal regression model
with pm.Model() as model:
    # Priors
    alpha = pm.Normal('alpha', mu=0, sigma=10)
    beta = pm.Normal('beta', mu=0, sigma=10)
    sigma = pm.HalfNormal('sigma', sigma=1)
    
    # Likelihood
    mu = alpha + beta * x_data
    y_obs = pm.Normal('y_obs', mu=mu, sigma=sigma, observed=y_data)
    
    # Sample with NUTS (4 chains, 2000 draws per chain)
    trace = pm.sample(2000, chains=4, return_inferencedata=True)

# Check convergence
pm.plot_trace(trace)
print(pm.summary(trace))

JAGS, WinBUGS, and OpenBUGS

JAGS (Just Another Gibbs Sampler), written by Martyn Plummer at IARC, is the modern successor to WinBUGS and OpenBUGS. It uses a BUGS-like modeling language and automatically constructs a Gibbs sampler (with Metropolis fallback) for any specified model. JAGS remains widely used in academic statistics and biostatistics, particularly for multilevel and hierarchical models. It’s called from R using the rjags or R2jags packages and from Python using pyjags. Generalized linear models are routinely fit with JAGS in Bayesian biostatistics courses.

WinBUGS (Bayesian inference Using Gibbs Sampling), developed at the MRC Biostatistics Unit in Cambridge, UK, pioneered accessible Bayesian computation in the early 1990s and directly enabled the Bayesian revolution in applied statistics. Its Windows-based interface lowered the barrier to MCMC considerably. OpenBUGS is the open-source continuation. Both are now largely superseded by Stan and PyMC in new projects, but remain important historically and are still used in some clinical trial and epidemiology settings in the UK and US. Clinical and healthcare research was among the earliest adopters of WinBUGS for Bayesian analysis of medical data.

Other Notable Implementations

Nimble is an R package developed at UC Berkeley that provides flexible MCMC algorithms with user-customizable samplers. It uses the BUGS modeling language but allows fine-grained control over which sampler is used for each parameter — useful for tailoring MCMC to specific model structures. TensorFlow Probability and NumPyro (JAX-based) bring MCMC to the deep learning ecosystem, enabling GPU-accelerated Bayesian neural networks and large-scale probabilistic models. Regularization in machine learning from a Bayesian perspective — where regularization equals a prior on weights — is implemented via these TensorFlow Probability and NumPyro frameworks.

Where MCMC Is Used: Applications Across Science and Industry

Markov Chain Monte Carlo has penetrated virtually every domain that requires inference from complex probabilistic models. Its spread from computational physics and statistics into biology, machine learning, economics, and social science reflects both its fundamental utility and the dramatic improvement in computing power since the 1990s. This section surveys the most important application domains — not just to illustrate range, but because understanding applications deepens intuition about why specific MCMC challenges (high dimensions, complex correlations, multimodality) arise and how they’re addressed.

Genomics and Computational Biology

Phylogenetics — the inference of evolutionary relationships from DNA sequences — is arguably one of the most important MCMC applications in biology. The software MrBayes, developed by John Huelsenbeck and Fredrik Ronquist, uses MCMC to sample from the joint posterior distribution over phylogenetic tree topologies, branch lengths, and substitution model parameters. The parameter space is a combination of discrete (tree topology) and continuous (branch lengths) components, making this a non-standard MCMC problem solved using reversible-jump MCMC. The journal Bioinformatics (Huelsenbeck & Ronquist, 2001) published the definitive MrBayes reference. Phylogenetic tree analysis in computational biology courses increasingly requires MCMC-based methods.

Genome-wide association studies (GWAS) and Bayesian fine-mapping of causal genetic variants also rely heavily on MCMC. With hundreds of thousands of genetic variants and complex linkage disequilibrium patterns, the posterior over causal variant configurations is high-dimensional and analytically intractable — exactly the setting where MCMC is indispensable. Researchers at the Wellcome Sanger Institute in Cambridge and the Broad Institute at MIT/Harvard use MCMC extensively for these analyses.

Epidemiology and Public Health

The COVID-19 pandemic brought Bayesian MCMC to global public attention through infectious disease modeling. The SIR (Susceptible-Infected-Recovered) and its variants use MCMC to infer transmission rates, reproduction numbers (Rₜ), and case ascertainment rates from noisy surveillance data. The Imperial College London** MRC Centre for Global Infectious Disease Analysis, led by researchers including Neil Ferguson, used Bayesian MCMC models extensively in pandemic response planning. The R package EpiEstim and the Python package PyMC-epidemiology provide MCMC-based R₀ estimation tools used by public health agencies worldwide. Causal inference methods in epidemiology combine with MCMC to address confounding and selection bias in observational disease studies.

Machine Learning and Probabilistic AI

Bayesian approaches to machine learning use MCMC to provide uncertainty quantification alongside predictions. Bayesian neural networks — where weights are treated as distributions rather than point estimates — use MCMC (typically HMC via NumPyro or TensorFlow Probability) to infer weight posteriors. This produces predictions with principled uncertainty bounds rather than just point estimates. Gaussian processes, widely used for spatial modeling, time series, and non-parametric regression, use MCMC to infer hyperparameters when conjugate inference isn’t available. Machine learning fundamentals increasingly include Bayesian and probabilistic treatments where MCMC is the computational backbone.

Latent Dirichlet Allocation (LDA) — a topic model for text analysis — was originally implemented using collapsed Gibbs sampling, exploiting the Dirichlet-multinomial conjugacy structure. LDA was developed by David Blei, Andrew Ng, and Michael Jordan at UC Berkeley and published in the Journal of Machine Learning Research (2003). The Gibbs sampler at its heart is a direct application of the coordinate-wise MCMC approach described in Section 3. Classification analysis methods like SVMs and decision trees stand in contrast to these Bayesian approaches — MCMC enables the uncertainty quantification that point-estimate classifiers cannot provide.

Economics and Social Science

Bayesian Vector Autoregression (BVAR) models are the dominant forecasting tool at central banks including the Federal Reserve, the Bank of England, and the European Central Bank. BVAR models use MCMC (typically Gibbs sampling) to infer joint posteriors over VAR coefficients and error covariances. The Minnesota prior — developed at the Federal Reserve Bank of Minneapolis — provides a hierarchical shrinkage prior that regularizes the high-dimensional parameter space. Without MCMC, fitting Bayesian macro-econometric models at the scale required for serious policy forecasting would be computationally impossible. Time series analysis with ARIMA methods are the classical alternative, but BVAR with MCMC provides richer uncertainty quantification.

Clinical Trials and Drug Development

Bayesian adaptive clinical trial designs increasingly use MCMC to update posterior probabilities of treatment effectiveness as data accumulates. The FDA (Food and Drug Administration) in the United States and the EMA (European Medicines Agency) both now accept Bayesian statistical analyses in regulatory submissions, with MCMC as the standard computation method. Berry Consultants, a statistical consulting firm specializing in adaptive trial design, deploys MCMC models for complex multi-arm trial analyses. Confidence intervals in classical statistics are replaced by posterior credible intervals in these Bayesian trial analyses — a conceptually important distinction with direct regulatory implications.

Advanced MCMC Methods: Beyond the Basics

The classical MCMC methods — Metropolis-Hastings and Gibbs sampling — are powerful but have limitations in specific problem structures. A rich ecosystem of advanced techniques has developed to address these limitations. Students and practitioners working on cutting-edge Bayesian inference should be aware of these methods, even if they rely on software implementations rather than coding them from scratch. Principal component analysis and reparametrization techniques connect directly to advanced MCMC methods, as both exploit structure in the parameter space for efficiency.

Slice Sampling

Slice sampling, developed by Radford Neal at the University of Toronto and published in the Annals of Statistics (2003), is a method for sampling from a one-dimensional distribution by “slicing” under the density curve. It adapts automatically to the local scale of the distribution without requiring explicit proposal tuning — a significant advantage over random walk Metropolis. Slice sampling is used within Gibbs samplers for full conditionals that don’t have closed forms, and it appears in NUTS as a component of the tree-building algorithm. It’s particularly useful for unimodal distributions with varying scale.

Parallel Tempering (Replica Exchange MCMC)

Parallel tempering (also called replica exchange MCMC) addresses multimodal posteriors by running multiple chains simultaneously at different “temperatures.” Higher-temperature chains have flattened posteriors and can easily traverse low-probability barriers between modes. Periodically, chains at adjacent temperatures propose to swap states — allowing the cold (target-temperature) chain to jump between modes via the high-temperature chains. Developed by Geyer and others in the early 1990s, parallel tempering is standard in computational chemistry for conformational sampling of molecular systems and in phylogenetics. It’s computationally expensive (k chains run simultaneously) but can solve multimodal problems where standard MCMC fails completely. Sampling distributions theory extends naturally to the multi-temperature chain framework of parallel tempering.

Reversible Jump MCMC (RJMCMC)

Reversible jump MCMC, developed by Peter Green at the University of Bristol and published in Biometrika (1995), extends MCMC to spaces of varying dimension. This is essential for problems where the number of parameters itself is unknown — for example, mixture models where the number of components, the order of a time series model, or the topology of a phylogenetic tree are all part of the posterior. RJMCMC allows jumps between parameter spaces of different dimensionality while maintaining detailed balance. It’s used in MrBayes for tree topology sampling and in Bayesian variable selection where the model dimension varies. Model selection with AIC and BIC is the frequentist alternative to the Bayesian model averaging that RJMCMC enables.

Sequential Monte Carlo (SMC)

Sequential Monte Carlo (SMC) — also called particle filters — is related to MCMC but works differently. Instead of one chain, SMC maintains a population of particles (samples) that are reweighted and resampled as data is processed sequentially. SMC is particularly suited to online inference (updating estimates as new data arrives) and state-space models (time series with latent states). Arnaud Doucet and colleagues at Oxford University have been central to SMC’s development, with the foundational textbook on Sequential Monte Carlo methods (2001) being the standard reference. PyMC implements SMC as an alternative to NUTS, particularly useful for tempered likelihoods and multimodal problems. Time series and state-space models are the primary setting where SMC outperforms standard MCMC.

Variational Inference: When MCMC Is Too Slow

Variational inference (VI) is not an MCMC method, but it’s important to understand as an alternative. Rather than sampling from the posterior, VI approximates it with a simpler parametric distribution, minimizing the KL divergence between the approximation and the true posterior. This is much faster than MCMC but less accurate — VI tends to underestimate posterior variance and can miss multimodality entirely. Automatic Differentiation Variational Inference (ADVI), developed by Kucukelbir and colleagues at Columbia, is implemented in Stan and PyMC and provides fast approximate posterior inference. Kucukelbir et al. (2017) in the Journal of Machine Learning Research is the key reference. The choice between MCMC and VI involves a speed-accuracy tradeoff that depends on the application. Overfitting and underfitting in machine learning have direct analogs in VI, where the choice of approximating family determines the bias-variance tradeoff.

Key People, Organizations, and Institutions in MCMC Development

The history and current state of MCMC has been shaped by a small number of pivotal figures and institutions. Knowing who they are and what makes each contribution unique is essential for academic writing on this topic — and it’s the kind of contextual depth that separates a sophisticated statistics essay from a surface-level summary. Writing research papers in statistics requires exactly this kind of entity-aware, citation-grounded exposition.

Nicholas Metropolis — Los Alamos National Laboratory

Nicholas Metropolis (1915–1999) was a Greek-American mathematician and physicist at Los Alamos National Laboratory, New Mexico. His 1953 paper, co-authored with Arianna and Marshall Rosenbluth, Augusta and Edward Teller, introduced what is now called the Metropolis algorithm. What makes Metropolis uniquely significant is context: the algorithm was born in the crucible of the Manhattan Project’s computational needs, where the challenge was computing equilibrium properties of thermodynamic systems too complex for analytic solutions. Metropolis built and programmed MANIAC I, one of the first electronic computers, specifically to run these simulations. His contribution is unique not just in the idea but in the practical synthesis of the computational hardware, the algorithmic invention, and the physical problem — all at once.

W.K. Hastings — University of Toronto

W. Keith Hastings (1930–2016) was a statistician at the University of Toronto who in 1970 published the crucial generalization of the Metropolis algorithm now known as the Metropolis-Hastings algorithm. What makes Hastings’ contribution uniquely significant is that he recognized the algorithm’s statistical generality — that the Metropolis algorithm was not just a physics tool but a general-purpose method for sampling from any distribution defined up to a normalizing constant. By allowing asymmetric proposal distributions and deriving the correct acceptance probability for them, Hastings transformed a specialized physics algorithm into the universal sampling method that powers modern Bayesian computation. Hastings’ 1970 paper in Biometrika is one of the most influential in statistical computing.

Andrew Gelman — Columbia University

Andrew Gelman is a Professor of Statistics and Political Science at Columbia University in New York City, and one of the most influential statisticians of the modern era. His contributions to MCMC include co-developing the Gelman-Rubin R-hat convergence diagnostic (with Donald Rubin, 1992), leading the development of the Stan probabilistic programming language and the No-U-Turn Sampler (with Matt Hoffman), and writing the definitive textbook on Bayesian Data Analysis (now in its third edition, co-authored with Carlin, Stern, Dunson, Vehtari, and Rubin). What makes Gelman uniquely significant is the breadth of his impact: he bridges theoretical statistics, practical computation, and applications in social science and epidemiology, and his public writing (including the Statistical Modeling, Causal Inference blog) has shaped how a generation of practitioners think about Bayesian inference and MCMC.

Radford Neal — University of Toronto

Radford Neal is a Professor at the University of Toronto whose contributions to MCMC are foundational. His 1996 PhD thesis and subsequent book “Bayesian Learning for Neural Networks” first applied HMC to neural network inference. His development of slice sampling (2003) provided an automatically-scaling alternative to M-H. His writing on MCMC theory, particularly on the importance of mixing and the geometry of the typical set, has been highly influential on algorithm design. What makes Neal uniquely significant is that he has repeatedly identified the deep geometric reasons why naive MCMC fails and invented algorithms — HMC, slice sampling — that directly address those geometric problems. The current HMC methods in Stan and PyMC trace directly to Neal’s theoretical and algorithmic innovations.

The Stan Development Team — Columbia University and Beyond

The Stan Development Team, originally centered at Columbia University and now an international open-source community, has built the most influential MCMC software ecosystem in statistics. Led by Andrew Gelman, Bob Carpenter, Matt Hoffman, and Michael Betancourt (among many others), Stan uniquely combines the NUTS algorithm, automatic differentiation, a rich modeling language, and comprehensive diagnostics into a single platform. What makes Stan uniquely significant is the rigorously Bayesian workflow it encodes — not just a sampler but a system that actively helps users understand and fix MCMC failures through divergence diagnostics, energy plots, and pair plots that reveal problematic posterior geometry.

Writing MCMC Assignments: Strategies for Statistics Students

MCMC appears in statistics and machine learning curricula in multiple forms — as theoretical questions about Markov chains and convergence, as applied exercises implementing algorithms in R or Python, and as analytical assignments interpreting output from Stan or PyMC models. Each type demands a different skill set. Statistics assignment help for MCMC topics requires integrating mathematical theory, computational skill, and clear written exposition — all three simultaneously. Research paper writing skills are especially relevant for MCMC assignments that require you to interpret posterior outputs and draw substantive conclusions.

Theoretical MCMC Questions: What Professors Want

Theoretical MCMC questions typically ask you to prove or demonstrate properties of specific algorithms. Common question types: prove that the Metropolis-Hastings algorithm satisfies detailed balance; derive the acceptance probability for a specific proposal distribution; show that a given Gibbs sampler is irreducible and aperiodic; compute the stationary distribution of a small discrete Markov chain. For these questions, the key is precision. Write your arguments formally using the language of probability — probability densities, transition kernels, conditional distributions. Don’t confuse “the chain will converge” with a proof that it has the right stationary distribution. Hypothesis testing theory and probability theory are the formal prerequisites for these proofs.

Computational MCMC Assignments: Implementation Quality

When asked to implement MCMC in code, professors evaluate three things: correctness (does the chain actually sample from the right distribution?), quality (is the implementation efficient, readable, and well-documented?), and diagnostic practice (do you actually check convergence?). A common mistake: implementing the algorithm and reporting results without any convergence checks. This is automatically penalized — it suggests you don’t understand that MCMC results are invalid without convergence verification. Always include trace plots, R-hat values, and ESS in your reported results. Choosing the right statistical test parallels choosing the right MCMC algorithm — the reasoning process is the same: understand the problem structure before selecting the method.

Interpreting MCMC Output for Applied Assignments

Applied Bayesian assignments — “fit this model to these data and interpret the results” — require you to go beyond just running the sampler. You need to: describe the model structure (priors, likelihood, parameters); verify convergence (R-hat, ESS, trace plots); report posterior summaries (means, credible intervals, tail probabilities); interpret those summaries substantively; and evaluate model fit using posterior predictive checks. Assignment rubric analysis for Bayesian statistics courses almost always awards marks specifically for posterior predictive checking and convergence diagnostics — sections that students frequently omit. Literature review writing for MCMC assignments should cite the original algorithm papers (Metropolis 1953, Hastings 1970, Gelfand & Smith 1990) alongside the software documentation you used.

Common MCMC Assignment Mistakes

The most common errors in MCMC assignments follow predictable patterns. Not running multiple chains — running a single chain makes it impossible to compute R-hat and easy to miss multimodal posteriors. Always run at least four chains in parallel. Reporting only the posterior mean — the whole point of MCMC is to characterize the full posterior distribution; reporting only the mean discards most of the information. Not tuning the proposal variance in M-H implementations — acceptance rates far from the 23% optimum (for Gaussian RW) signal poor exploration. Interpreting credible intervals as confidence intervals — they have different interpretations; a 95% credible interval means the parameter is in that range with 95% posterior probability, not that 95% of such intervals would contain the true value. Common assignment mistakes in statistics are often about precision of language — and the credible vs. confidence interval distinction is a prime example.

MCMC Glossary: Essential Terms, LSI Keywords, and Related Concepts

Mastery of MCMC in academic writing, exams, and professional practice requires command of the field’s precise vocabulary. The following terms are those most frequently tested in statistics courses, appearing in rubrics, exam questions, and peer-reviewed literature. They are also the natural language processing (NLP) and latent semantic indexing (LSI) terms that define the conceptual landscape around this topic. Probability distributions and random variables are the foundational mathematical objects from which all MCMC concepts are built.

Core MCMC Vocabulary

Stationary distribution — the target distribution that the Markov chain is designed to sample from; once reached, the chain’s distribution doesn’t change. Ergodicity — the property ensuring a chain will eventually visit all regions of the target distribution, regardless of starting point. Detailed balance — the reversibility condition π(x)K(x, x’) = π(x’)K(x’, x) that guarantees a chain’s stationary distribution is π. Mixing time — how long the chain takes to effectively forget its starting point. Proposal distribution — the distribution from which candidate moves are drawn in Metropolis-Hastings. Acceptance probability — the probability of accepting a proposed move, ensuring the chain samples proportionally from the target. Burn-in — initial chain iterations before convergence, to be discarded. Thinning — retaining only every kth sample to reduce storage; does not increase ESS.

Trace plot — visualization of parameter values across iterations; well-mixed chains look like “fuzzy caterpillars.” R-hat (Gelman-Rubin statistic) — convergence diagnostic comparing within- and between-chain variance; values near 1.0 indicate convergence. Effective sample size (ESS) — the equivalent number of independent samples, accounting for autocorrelation. Autocorrelation — correlation between successive samples; high autocorrelation means low ESS. Divergent transition — HMC-specific failure indicating numerical instability in the leapfrog integrator, signaling posterior geometry problems. Warmup — Stan’s adaptive burn-in phase where proposal parameters are tuned automatically. Posterior predictive check — comparing simulated data from the fitted model to observed data for model evaluation. Residual analysis in regression serves an analogous model-checking function in frequentist frameworks.

Related Statistical Concepts and LSI Terms

Bayesian inference — the inferential framework using Bayes’ theorem to update prior beliefs with data; MCMC is its primary computational tool. Posterior distribution — the target of Bayesian inference; prior × likelihood, normalized. Likelihood function — the probability of observed data as a function of model parameters. Prior distribution — encodes parameter beliefs before data. Conjugate prior — a prior that yields a posterior in the same distributional family; enables closed-form Gibbs sampling. Hierarchical model — a model where parameters are drawn from group-level distributions, themselves drawn from hyperpriors; the natural setting for MCMC. Marginal likelihood — the normalizing constant in Bayes’ theorem; MCMC sidesteps its computation. Credible interval — Bayesian uncertainty interval: P(a < θ < b | data) = 0.95. Highest density interval (HDI) — the narrowest credible interval containing a specified probability mass.

Related methods and concepts that appear in MCMC literature: importance sampling (reweighted direct sampling), rejection sampling (proposal-accept direct sampling), Langevin dynamics (MCMC using gradient information without full HMC), particle filter (sequential MC for time series), variational Bayes (approximate posterior inference), expectation-maximization (EM) (frequentist analog for latent variable models), annealing (temperature-based MCMC for multimodal targets), tempering (parallel chains at different temperatures). Missing data imputation methods using multiple imputation are essentially MCMC-based, as the imputations represent samples from the posterior predictive distribution of missing values.

For students writing literature reviews or extensive analyses involving MCMC, the most important journals are: Journal of the American Statistical Association (JASA), Annals of Statistics, Biometrika, Journal of Applied Bayesian Analysis, Statistical Science, Journal of Machine Learning Research (JMLR), and Journal of Chemical Physics (for the historical computational physics origins). Writing an exemplary literature review for an MCMC topic requires navigating sources from statistics, computer science, physics, and domain-specific journals — a genuinely cross-disciplinary challenge.

Frequently Asked Questions: Markov Chain Monte Carlo (MCMC)