Introduction

Imagine you want to know the average height of a mountain range — but you can only take measurements at random spots, and the range is so vast and jagged that you cannot compute the average analytically. You wander around, spending more time in the tall regions and less in the valleys, and eventually your sample reflects the shape of the terrain. That wandering is, in essence, Markov Chain Monte Carlo.

MCMC is a family of algorithms for sampling from a probability distribution when direct calculation is impossible. The workhorse is the Metropolis algorithm, invented in 1953 by Nicholas Metropolis, Arianna Rosenbluth, Marshall Rosenbluth, Augusta Teller, and Edward Teller while simulating the thermodynamics of hard-sphere molecules. Decades later it was recognized as one of the ten most important algorithms of the twentieth century.

The core insight is disarmingly simple: you do not need to know everything about a distribution to sample from it. You only need to compare relative probabilities — is the proposed next state more or less likely than the current one? That ratio is almost always easy to compute, even when the normalizing constant (the total area under the curve) is completely intractable.

This single trick unlocks Bayesian inference, statistical physics, molecular dynamics, and modern machine learning. It is solved algorithmically — the Metropolis–Hastings algorithm is proven to converge to the correct distribution in the limit — but the practical challenge of mixing time (how long the walk must run before its samples are trustworthy) remains an active research frontier.

Watch the Walk

The target distribution below has two peaks — a mixture of two Gaussians. A naive approach would require computing the normalizing constant; the Metropolis algorithm ignores it entirely.

Press Run to start the random walk. At each step a candidate position is proposed; if it is more likely than the current position it is always accepted; if less likely it is accepted with probability proportional to how much less likely it is. The histogram below the canvas shows where the walk has spent its time.

<div class="controls">
  <label>Step size <span id="stepVal">1.0</span>
    <input type="range" id="stepSize" min="0.1" max="4" step="0.1" value="1.0">
  </label>
  <button id="runBtn" type="button">Run</button>
  <button id="stepBtn" type="button">Step ×10</button>
  <button id="resetBtn" type="button" class="ghost">Reset</button>
  <span class="stat">Samples: <b id="nSamples">0</b> &nbsp; Accept: <b id="acceptRate">—</b></span>
</div>
<canvas id="canvas" width="560" height="200"></canvas>
<div id="histWrap"><canvas id="hist" width="560" height="120"></canvas></div>
<p class="hint">The <b>curve</b> is the target (bimodal Gaussian). The <b>dot</b> is the current walker. The <b>histogram</b> accumulates visited positions — it should match the curve.</p>

* { box-sizing: border-box; }
body { font-family: system-ui, sans-serif; color: #222; margin: 0; background: #fff; }
.controls { display: flex; flex-wrap: wrap; align-items: center; gap: .5rem .8rem; margin-bottom: .5rem; }
label { font-size: .85rem; display: flex; align-items: center; gap: .4rem; }
input[type=range] { width: 90px; }
button { font: 600 13px system-ui; padding: .38rem .8rem; border: 1px solid #1d3557;
         background: #1d3557; color: #fff; border-radius: 8px; cursor: pointer; }
button.ghost { background: #fff; color: #1d3557; }
.stat { font-size: .83rem; color: #444; }
canvas { display: block; border-radius: 6px; background: #f4f7fa; width: 100%; height: auto; }
#histWrap { margin-top: 4px; }
.hint { font-size: .8rem; color: #555; margin: .4rem 0 0; line-height: 1.4; }

// Metropolis algorithm on a bimodal Gaussian target:
//   p(x) ∝ 0.4·N(x; -2, 0.8²) + 0.6·N(x; 2, 1²)
// We never need the normalizing constant — only ratios p(x')/p(x).

const BINS = 56, X_MIN = -6, X_MAX = 6;
const W = 560, H = 200, HH = 120;

// --- Target distribution (unnormalized) ---
function gauss(x, mu, sigma) {
  return Math.exp(-0.5 * ((x - mu) / sigma) ** 2) / sigma;
}
function target(x) {
  return 0.4 * gauss(x, -2, 0.8) + 0.6 * gauss(x, 2, 1.0);
}

// Precompute the target curve for drawing
const CURVE_SAMPLES = 300;
const curveX = [], curveY = [];
let curveMax = 0;
for (let i = 0; i <= CURVE_SAMPLES; i++) {
  const x = X_MIN + (X_MAX - X_MIN) * i / CURVE_SAMPLES;
  const y = target(x);
  curveX.push(x); curveY.push(y);
  if (y > curveMax) curveMax = y;
}

// --- State ---
let current = 0;          // current position of the chain
let counts = new Array(BINS).fill(0);
let nSamples = 0;
let nAccepted = 0;
let running = false;
let rafId = null;
let stepSize = 1.0;

// --- DOM ---
const canvas  = document.getElementById('canvas');
const ctx     = canvas.getContext('2d');
const histCvs = document.getElementById('hist');
const hctx    = histCvs.getContext('2d');
const runBtn  = document.getElementById('runBtn');
const stepBtn = document.getElementById('stepBtn');
const resetBtn= document.getElementById('resetBtn');
const stepSlider = document.getElementById('stepSize');
const stepVal = document.getElementById('stepVal');
const nEl     = document.getElementById('nSamples');
const accEl   = document.getElementById('acceptRate');

stepSlider.addEventListener('input', () => {
  stepSize = parseFloat(stepSlider.value);
  stepVal.textContent = stepSize.toFixed(1);
});

// --- Coordinate helpers ---
function xToPixel(x)  { return (x - X_MIN) / (X_MAX - X_MIN) * W; }
function yToPixel(y, maxY, canvasH) { return canvasH - y / maxY * (canvasH - 12) - 6; }

// --- Metropolis step ---
function metropolisStep() {
  const proposal = current + (Math.random() * 2 - 1) * stepSize;
  const ratio = target(proposal) / target(current);
  if (Math.random() < ratio) {
    current = proposal;
    nAccepted++;
  }
  // bin the current position
  const bin = Math.floor((current - X_MIN) / (X_MAX - X_MIN) * BINS);
  if (bin >= 0 && bin < BINS) counts[bin]++;
  nSamples++;
}

// --- Draw target curve + walker ---
function drawCanvas() {
  ctx.clearRect(0, 0, W, H);

  // filled area under curve
  ctx.beginPath();
  ctx.moveTo(xToPixel(curveX[0]), yToPixel(curveY[0], curveMax, H));
  for (let i = 1; i <= CURVE_SAMPLES; i++) {
    ctx.lineTo(xToPixel(curveX[i]), yToPixel(curveY[i], curveMax, H));
  }
  ctx.lineTo(xToPixel(curveX[CURVE_SAMPLES]), H);
  ctx.lineTo(xToPixel(curveX[0]), H);
  ctx.closePath();
  ctx.fillStyle = 'rgba(29,53,87,0.10)';
  ctx.fill();

  // curve stroke
  ctx.beginPath();
  ctx.moveTo(xToPixel(curveX[0]), yToPixel(curveY[0], curveMax, H));
  for (let i = 1; i <= CURVE_SAMPLES; i++) {
    ctx.lineTo(xToPixel(curveX[i]), yToPixel(curveY[i], curveMax, H));
  }
  ctx.strokeStyle = '#1d3557';
  ctx.lineWidth = 2;
  ctx.stroke();

  // walker dot
  const px = xToPixel(current);
  const py = yToPixel(target(current), curveMax, H);
  ctx.beginPath();
  ctx.arc(px, py, 7, 0, Math.PI * 2);
  ctx.fillStyle = '#e63946';
  ctx.fill();

  // x-axis labels
  ctx.fillStyle = '#888';
  ctx.font = '11px system-ui';
  ctx.textAlign = 'center';
  for (let x = X_MIN; x <= X_MAX; x += 2) {
    ctx.fillText(x, xToPixel(x), H - 2);
  }
}

// --- Draw histogram ---
function drawHist() {
  hctx.clearRect(0, 0, W, HH);
  const maxCount = Math.max(...counts, 1);
  const binW = W / BINS;
  for (let i = 0; i < BINS; i++) {
    const barH = counts[i] / maxCount * (HH - 16);
    hctx.fillStyle = 'rgba(29,53,87,0.55)';
    hctx.fillRect(i * binW + 1, HH - 14 - barH, binW - 2, barH);
  }
  // overlay target curve on histogram
  hctx.beginPath();
  for (let i = 0; i <= CURVE_SAMPLES; i++) {
    const px = xToPixel(curveX[i]);
    const py = HH - 14 - curveY[i] / curveMax * (HH - 16);
    i === 0 ? hctx.moveTo(px, py) : hctx.lineTo(px, py);
  }
  hctx.strokeStyle = '#e63946';
  hctx.lineWidth = 2;
  hctx.stroke();

  // axis
  hctx.fillStyle = '#999';
  hctx.font = '10px system-ui';
  hctx.textAlign = 'center';
  for (let x = X_MIN; x <= X_MAX; x += 2) {
    hctx.fillText(x, xToPixel(x), HH - 2);
  }
}

function updateStats() {
  nEl.textContent = nSamples.toLocaleString();
  accEl.textContent = nSamples > 0
    ? (nAccepted / nSamples * 100).toFixed(1) + '%'
    : '—';
}

function render() {
  drawCanvas();
  drawHist();
  updateStats();
}

// --- Animation loop ---
const STEPS_PER_FRAME = 30;
function loop() {
  for (let i = 0; i < STEPS_PER_FRAME; i++) metropolisStep();
  render();
  if (running) rafId = requestAnimationFrame(loop);
}

runBtn.addEventListener('click', () => {
  running = !running;
  runBtn.textContent = running ? 'Pause' : 'Run';
  if (running) rafId = requestAnimationFrame(loop);
  else cancelAnimationFrame(rafId);
});

stepBtn.addEventListener('click', () => {
  for (let i = 0; i < 10; i++) metropolisStep();
  render();
});

resetBtn.addEventListener('click', () => {
  running = false;
  runBtn.textContent = 'Run';
  cancelAnimationFrame(rafId);
  current = 0;
  counts = new Array(BINS).fill(0);
  nSamples = 0;
  nAccepted = 0;
  render();
});

render();
Result
Console

Notice how the walk occasionally tunnels from one peak to the other — hopping through the low-probability valley. The histogram converges toward the true target (shown as a curve). Too small a step size keeps the walk exploring only locally; too large a step size causes most proposals to overshoot low-probability regions and get rejected. Adjusting the step size slider reveals this tradeoff in real time.

The Real Complexity

The Metropolis algorithm is proven correct: run it long enough and the fraction of time spent near any point converges to that point's probability under the target distribution. This follows from the theory of ergodic Markov chains — chains that eventually reach every state and do not get trapped in cycles.

But "long enough" hides the difficulty:

  • Mixing time is the number of steps needed before the chain's distribution is close (in total variation distance) to the stationary distribution. For well-behaved, unimodal distributions this is polynomial in the dimension. For multimodal distributions separated by low-probability barriers, mixing can take exponentially long — the chain is stuck in one peak for an astronomically long time before accidentally crossing the valley.
  • Spectral gap: the speed of mixing is governed by the second-largest eigenvalue of the chain's transition matrix. A small spectral gap means slow mixing. Computing this gap is hard in general.
  • Dimension curse: in high dimensions, typical random-walk proposals land in low-probability regions almost always, so acceptance rates collapse. Algorithms like Hamiltonian Monte Carlo (HMC) use gradient information to make smarter proposals, keeping acceptance rates healthy even in hundreds of dimensions. HMC underlies modern probabilistic programming systems like Stan and PyMC.
  • No free lunch: there is no universal fast sampler. Deciding whether a given Markov chain mixes in polynomial time is related to questions about counting complexity (the #P class) and remains open in general.

The status is therefore nuanced: Metropolis–Hastings is solved as an algorithm (it converges), but efficient sampling from arbitrary distributions is not — it can require exponential time in the worst case, and proving that a specific chain mixes fast is often a hard open problem.

Where It Matters

MCMC's ability to sample from any distribution whose unnormalized density can be evaluated makes it indispensable across science and engineering:

  • Bayesian inference: posterior distributions in statistical models almost never have closed forms. MCMC is the standard engine for approximate Bayesian computation — from simple regressions to hierarchical models with thousands of parameters.
  • Molecular dynamics and protein folding: simulating how atoms explore energy landscapes to find stable configurations. The original Metropolis paper was literally about molecules.
  • Ising model and statistical physics: simulating phase transitions in magnets and lattice models — exactly the setting where the algorithm was born.
  • Image reconstruction: in medical imaging (MRI, CT), reconstructing an image from noisy measurements is a Bayesian problem; MCMC explores the space of plausible images.
  • Probabilistic machine learning: variational inference and Monte Carlo gradient estimators underpin modern deep learning. Even the training of diffusion models (the backbone of generative image AI) is intimately connected to Langevin dynamics — a continuous-time cousin of Metropolis.
  • Phylogenetics: inferring evolutionary trees from DNA sequence data is done almost exclusively with MCMC (the MrBayes and BEAST programs).

In each case the pattern is the same: the space of possibilities is too vast to enumerate, but relative probabilities are easy to compute. The random walk does the rest.

Conclusion

The Metropolis algorithm is one of the great algorithmic ideas of the twentieth century: a random walk that needs only relative probabilities, yet provably converges to sample from any target distribution. It turned intractable integrals into tractable simulations overnight, and the applications have never stopped growing.

The catch is mixing time — the walk must wander long enough to explore the full distribution, and in the worst case that takes exponentially long. Modern variants like Hamiltonian Monte Carlo and Langevin dynamics add gradient information to speed the exploration, but the fundamental tension between correctness and efficiency remains unsolved in general.

If you have encountered Bayesian inference or wondered how generative AI models explore their parameter spaces, you have brushed against MCMC. The random walk is everywhere — and it still has secrets left to give up.

Share this article

Pick a channel — or use your device's native share sheet.

Comments

Loading comments...

https://www.kipuhub.com/en/article/mcmc-metropolis/Content licensed under CC BY-NC 4.0.