MCMC Metropolis-Hastings in 50 Lines of JavaScript

1. What Is MCMC and Why Do We Need It

Suppose you know a probability density up to a constant, p(x) = f(x)/Z, where f(x) is easy to evaluate but the normalising constant Z = ∫ f(x) dx is intractable (a common situation in Bayesian inference, where f is prior × likelihood and Z is the marginal likelihood). You cannot sample from p(x) directly, but you can still evaluate f(x) at any point.

Markov Chain Monte Carlo (MCMC) constructs a Markov chain whose stationary distribution equals p(x) — a random walk that, after enough steps, visits each region of the state space with exactly the right long-run frequency, without ever needing to know Z. Metropolis-Hastings (1953/1970) is the classical, general-purpose recipe for building such a chain for essentially any target density.

2. Target and Proposal Distributions

The target distribution π(x) ∝ f(x) is what we want samples from. The proposal distribution q(x'|x) generates a candidate next state x' given the current state x. The simplest and most common choice is a symmetric random-walk proposal:

Because this proposal is symmetric (q(x'|x) = q(x|x')), it simplifies the acceptance formula in the next step — this special case is called the Metropolis algorithm. Metropolis- Hastings generalises it to asymmetric proposals with an extra correction term.

3. The Metropolis-Hastings Acceptance Rule

At each step, propose a candidate x' from q(x'|x), then accept it with probability α, otherwise stay at x:

For a symmetric proposal, q(x|x') = q(x'|x) cancels, leaving the simpler Metropolis ratio:

Intuition: if the candidate has higher density than the current state (f(x') > f(x)), always move there. If it has lower density, move there only sometimes, with probability equal to the density ratio — this occasional acceptance of "worse" moves is exactly what lets the chain explore the full distribution instead of climbing straight to the single mode and getting stuck.

4. The 50-Line Implementation

Below, logTarget(x) returns the log of the unnormalised target density (working in log-space avoids underflow for tiny probabilities and turns the ratio into a subtraction). The example target is a mixture of two Gaussians — a bimodal distribution that would be awkward to sample directly with Math.random()-based inverse-CDF methods:

// Unnormalised log-density: mixture of N(-2, 1) and N(3, 0.7), weights 0.4/0.6
function logTarget(x) {
  const logNormal = (x, mu, sigma) =>
    -0.5 * ((x - mu) / sigma) ** 2 - Math.log(sigma);
  const a = Math.log(0.4) + logNormal(x, -2, 1);
  const b = Math.log(0.6) + logNormal(x, 3, 0.7);
  // log-sum-exp for numerical stability: log(e^a + e^b)
  const m = Math.max(a, b);
  return m + Math.log(Math.exp(a - m) + Math.exp(b - m));
}

// Standard Normal sample via Box-Muller transform
function gaussianStep(sigma) {
  const u1 = Math.random(), u2 = Math.random();
  return sigma * Math.sqrt(-2 * Math.log(u1)) * Math.cos(2 * Math.PI * u2);
}

// Metropolis-Hastings sampler (symmetric random-walk proposal)
function metropolisHastings(logTarget, {
  nSamples = 20000, sigma = 1.0, x0 = 0,
} = {}) {
  const samples = new Array(nSamples);
  let x = x0;
  let logFx = logTarget(x);
  let accepted = 0;

  for (let i = 0; i < nSamples; i++) {
    const xProp   = x + gaussianStep(sigma);
    const logFxp  = logTarget(xProp);
    const logAlpha = logFxp - logFx;               // log of f(x')/f(x)

    if (Math.log(Math.random()) < logAlpha) {   // accept if log(U) < log(alpha), alpha capped at 1
      x = xProp;
      logFx = logFxp;
      accepted++;
    }
    // else: reject, stay at current x (this repeat IS the sample)
    samples[i] = x;
  }
  return { samples, acceptanceRate: accepted / nSamples };
}

const { samples, acceptanceRate } = metropolisHastings(logTarget, { nSamples: 20000, sigma: 1.2 });
console.log('acceptance rate:', acceptanceRate.toFixed(3)); // aim for ~0.2-0.5

That is the entire algorithm: propose, evaluate, accept-or-reject, repeat. Everything else in this tutorial — burn-in, thinning, diagnostics, multivariate targets — is refinement around this ~30-line core.

5. Burn-In and Thinning

The chain starts at an arbitrary x0, which may sit in a low-probability region far from where the target's mass actually is. The initial samples are therefore not representative of π(x) and must be discarded — this is the burn-in period.

function applyBurnInAndThin(samples, { burnIn = 2000, thin = 5 } = {}) {
  return samples
    .slice(burnIn)               // discard the first `burnIn` samples
    .filter((_, i) => i % thin === 0); // keep every `thin`-th sample
}

Thinning — keeping only every k-th sample — reduces the strong autocorrelation between consecutive samples of a random-walk chain (each sample is close to its predecessor), at the cost of throwing away most of the raw samples. It does not increase the true information content per kept sample, but it dramatically simplifies downstream code that assumes roughly-independent samples (e.g. naive histogram density estimates or simple summary statistics).

6. Diagnostics: Trace Plots and Autocorrelation

Never trust an MCMC run blindly — always check that it actually converged and mixed well:

• Trace plot: plot sample value against iteration number. A well-mixed chain looks like noisy "fuzzy caterpillar" oscillation around the target's support; a poorly-mixed one shows long flat segments (the chain got stuck) or a slow drift (still burning in).
• Acceptance rate: too high (> 70%) usually means the proposal step σ is too small — the chain crawls and mixes slowly. Too low (< 10%) means σ is too large — most proposals are rejected. The classical target for random-walk Metropolis in 1D is roughly 20–50%.
• Autocorrelation function (ACF): measures how correlated sample i is with sample i+k. Slow-decaying ACF implies a large effective sample size penalty.

// Lag-k autocorrelation of the raw (unthinned) chain
function autocorrelation(samples, maxLag = 50) {
  const n = samples.length;
  const mean = samples.reduce((a,b) => a+b, 0) / n;
  const variance = samples.reduce((s,x) => s + (x-mean)**2, 0) / n;
  const acf = [];
  for (let k = 0; k <= maxLag; k++) {
    let cov = 0;
    for (let i = 0; i < n - k; i++) cov += (samples[i]-mean) * (samples[i+k]-mean);
    acf.push(cov / (n - k) / variance);
  }
  return acf; // acf[0] === 1 always; watch how fast it decays toward 0
}

// Effective sample size (rough estimate via sum of ACF)
function effectiveSampleSize(samples) {
  const acf = autocorrelation(samples, 100);
  let sum = 1;
  for (let k = 1; k < acf.length; k++) {
    if (acf[k] < 0.05) break;   // stop once autocorrelation is negligible
    sum += 2 * acf[k];
  }
  return Math.round(samples.length / sum);
}

7. Extending to Multiple Dimensions and Tuning the Step Size

The core algorithm generalises directly to a vector state x ∈ ℝᵈ — replace the scalar Gaussian step with a multivariate Normal step, and logTarget with a function of a vector:

function gaussianStepVec(dim, sigma) {
  return Array.from({ length: dim }, () => gaussianStep(sigma));
}

function metropolisHastingsND(logTarget, dim, { nSamples = 20000, sigma = 0.5 } = {}) {
  let x = new Array(dim).fill(0);
  let logFx = logTarget(x);
  const samples = [];
  for (let i = 0; i < nSamples; i++) {
    const step = gaussianStepVec(dim, sigma);
    const xProp = x.map((v, j) => v + step[j]);
    const logFxp = logTarget(xProp);
    if (Math.log(Math.random()) < logFxp - logFx) { x = xProp; logFx = logFxp; }
    samples.push([...x]);
  }
  return samples;
}

As dimension d grows, a naive isotropic random-walk step becomes increasingly inefficient — the acceptance rate collapses unless σ shrinks roughly as 1/√d (a classical result for optimal scaling of random-walk Metropolis). In high dimensions, more advanced samplers such as Hamiltonian Monte Carlo (HMC) or the No-U-Turn Sampler (NUTS) — used by Stan and PyMC — use gradient information to propose far more efficient moves. Metropolis-Hastings remains the essential building block that explains why all of them work.