Metropolis–Hastings Sampler Simulation

The algorithm

Let π(x) be the target density, known only up to a constant. From the current state x:

1. Propose y = x + σ_prop·Z with Z ~ N(0, 1).
2. Compute the acceptance probability α = min(1, π(y)·q(x | y) / [π(x)·q(y | x)]). This Gaussian random walk is symmetric, q(x | y) = q(y | x), so α = min(1, π(y)/π(x)) — Metropolis's original rule; Hastings added the q ratio for asymmetric proposals.
3. Draw u ~ U(0, 1). If u < α, move to y; otherwise stay at x and record x again.

Only the ratio π(y)/π(x) is ever needed, so the normalizing constant — often an intractable integral in Bayesian statistics — cancels. The rule makes the chain satisfy detailed balance, π(x)P(x → y) = π(y)P(y → x), which makes π its stationary distribution. The code works with log densities to avoid underflow in the tails. Staying put after a rejection is essential: dropping the repeated values would bias the sample toward low-density regions.

Step size: too small, too large, about right

The chain starts at x₀ = 6, far in the tail of the standard normal, so the first dozen or so iterations are a transient that says nothing about the target. The shaded burn-in period is discarded. After that, three regimes are worth trying:

• Tiny steps (σ_prop = 0.1). Almost every proposal is accepted, yet the chain crawls: successive samples are nearly identical, the trace looks like a slow random walk and the autocorrelation bars barely fall over 40 lags. High acceptance is not the goal.
• Huge steps (σ_prop = 20). Most proposals land far in the tails and are rejected, so the trace shows long flat stretches. The samples are again strongly correlated.
• In between (σ_prop ≈ 2.4). For a one-dimensional normal target this is the classic rule of thumb, giving an acceptance rate near 44% and the largest effective sample size of the three.

The bimodal target shows a harder failure. Its two modes are four units apart with almost no probability between them, so with σ_prop = 1 the chain crosses only a few dozen times in 10,000 iterations. The share of time spent in each mode — and with it the sample mean — rests on those few crossings and swings noticeably from seed to seed, while the acceptance rate looks perfectly healthy. Only the trace and the known mean (0.8) give it away. Steps between about 2.4 and 4 cross hundreds of times. The Student t target with 3 degrees of freedom has variance 3 but an infinite fourth moment, so the sample standard deviation converges slowly and usually undershoots √3 ≈ 1.732 in short runs.

Effective sample size

Correlated samples carry less information than independent ones. If ρ_k is the autocorrelation at lag k, the variance of the sample mean of n draws is approximately σ²·τ/n, with integrated autocorrelation time τ = 1 + 2·Σ ρ_k. The effective sample size is

The sum cannot run over all lags, because the far-lag estimates are pure noise. The lab uses Geyer's initial positive sequence: it adds pairs Γ_m = ρ_2m + ρ_2m+1 while they are positive, so τ = −1 + 2·Σ Γ_m, with the autocorrelations computed from the kept samples (divisor n) and lags capped at 1,000. If the cap is reached, the status line marks the ESS as an upper bound. For a chain that behaves like an AR(1) process with lag-1 correlation ρ, τ ≈ (1 + ρ)/(1 − ρ) — a useful mental check against the lag-1 readout. Monte Carlo error bars from MCMC output should use the ESS, not the raw number of iterations; see how accurate a Monte Carlo estimate is.

Example result

With the standard normal target, σ_prop = 1, 10,000 iterations, 1,000 burn-in and seed 42, the chain accepts 71.2% of proposals. The 9,000 kept samples have mean 0.037 and standard deviation 1.015 (target 0 and 1), but a lag-1 autocorrelation of 0.782 shrinks their worth to an effective sample size of about 1,070.

Implementation

The page runs simulation.js from the package below; the Python file is a line-by-line port with the same random-number generator, so a given seed produces the same result in both languages. The excerpts show the core of each.

// Random-walk Metropolis: propose y = x + σ_prop·Z and accept with probability
// min(1, π(y)/π(x)). The Gaussian proposal is symmetric, so the Hastings
// correction q(x|y)/q(y|x) equals 1. Only ratios of π are needed: no normalising constant.
export function mhStep(x, logpx, target, step, rng) {
  const y = x + step * rng.normal();
  const logpy = target.logp(y);
  const u = rng.uniform();
  const logA = logpy - logpx;
  if (logA >= 0 || u < Math.exp(logA)) return { x: y, logp: logpy, accepted: true };
  return { x, logp: logpx, accepted: false };
}

export function runChain(target, start, step, iterations, rng) {
  const chain = new Float64Array(iterations + 1);
  const accepted = new Uint8Array(iterations + 1);
  let s = { x: start, logp: target.logp(start), accepted: false };
  chain[0] = start;
  for (let i = 1; i <= iterations; i++) {
    s = mhStep(s.x, s.logp, target, step, rng);
    chain[i] = s.x;
    accepted[i] = s.accepted ? 1 : 0;
  }
  return { chain, accepted };
}

def run_chain(logp, start, step, iterations, rng):
    """Random-walk Metropolis: propose y = x + step*Z, accept with probability min(1, pi(y)/pi(x))."""
    x, lpx = start, logp(start)
    chain, accepted = [start], [0]
    for _ in range(iterations):
        y = x + step * rng.normal()
        lpy = logp(y)
        u = rng.uniform()
        log_a = lpy - lpx
        if log_a >= 0 or u < math.exp(log_a):  # symmetric proposal: no Hastings correction
            x, lpx = y, lpy
            accepted.append(1)
        else:
            accepted.append(0)
        chain.append(x)
    return chain, accepted

Run the full model from a terminal:

node simulation.js --target normal --step 1 --seed 42
python simulation.py --target normal --step 1 --seed 42

Download the model

Assumptions & limitations

• Random-walk proposal only. Independence samplers, Hamiltonian Monte Carlo and Gibbs sampling are other Metropolis–Hastings variants with very different behavior, especially in many dimensions.
• One dimension. The targets are one-dimensional so the histogram can be checked by eye; in d dimensions the good random-walk step shrinks roughly like 1/√d, and the ideal acceptance rate drops toward 23%.
• Burn-in is a fixed count. Discarding a fixed number of iterations, capped at half the run, is a heuristic. There is no general test that proves a chain has converged; multiple chains from dispersed starting points are the usual safeguard.
• ESS is an estimate. The initial-positive-sequence estimator is itself noisy, and for chains stuck in one mode it can report a comfortable ESS for the wrong distribution.
• Pseudo-random numbers. Proposals and acceptance draws come from a seeded generator; see pseudo-random numbers and seeds.

Related simulations

• OS-014Markov chainA process that remembers only its current state: stationary distribution, mixing and absorption, live.
• OS-002Monte Carlo integrationEstimate an integral as the average of f at random points, with an error bar you can compute.
• OS-038Ising modelSpins on a grid that like to align: a phase transition at T_c ≈ 2.269, simulated with Metropolis flips.
• OS-012Random walkCoin-flip steps on a line or grid: distance grows like √n, and in 1-D and 2-D every walker comes back.

Background: How Accurate Is a Monte Carlo Estimate? · Pseudo-Random Numbers and Seeds in Simulation

References

1. Metropolis, N., Rosenbluth, A. W., Rosenbluth, M. N., Teller, A. H. & Teller, E. (1953). Equation of state calculations by fast computing machines. Journal of Chemical Physics, 21(6), 1087–1092. doi.org/10.1063/1.1699114
2. Hastings, W. K. (1970). Monte Carlo sampling methods using Markov chains and their applications. Biometrika, 57(1), 97–109. doi.org/10.1093/biomet/57.1.97
3. Robert, C. P. & Casella, G. (2004). Monte Carlo Statistical Methods, 2nd ed. Springer.