Metropolis-Hastings Sampling a Gaussian
Draw samples from a target distribution using a random-walk proposal and the accept-reject rule, then check the sample statistics.
Problem
Metropolis-Hastings is a Markov chain Monte Carlo method that samples from a distribution known only up to a constant. It proposes a move, then accepts or rejects it with a probability that guarantees the chain's stationary distribution is the target. We sample a standard Gaussian to keep the check simple.
Algorithm
From the current point x, propose x' = x + normal(0, step). Accept with probability min(1, p(x')/p(x)); for a Gaussian target that ratio is exp(-(x'^2 - x^2)/2). Rejected moves keep the current point, which is what preserves the correct distribution.
import numpy as np
rng=np.random.default_rng(11)
def logp(x): return -0.5*x*x
x=0.0; step=1.5; samples=[]
for _ in range(50000):
xp=x+rng.normal(0,step)
if np.log(rng.random())<logp(xp)-logp(x):
x=xp
samples.append(x)
s=np.array(samples[1000:])
print('mean',round(s.mean(),3),'std',round(s.std(),3))
Result
After discarding an initial burn-in the sample mean is near 0 and the standard deviation near 1, matching the target Gaussian. The step size tunes efficiency: too small and the chain crawls, too large and most proposals are rejected. A common target acceptance rate is around 25 to 45 percent. Successive samples are correlated, so effective sample size is smaller than the raw count.
- Only the ratio of densities is needed, so the normalizing constant never has to be computed.
- Convergence requires discarding burn-in and checking mixing diagnostics before trusting estimates.
- MCMC underlies Bayesian parameter inference in Kronos uncertainty-quantification studies.