Understanding the Metropolis-Hastings Algorithm: A Step-by-Step Guide with Python Examples

The Metropolis-Hastings (MH) algorithm is a cornerstone of Markov Chain Monte Carlo (MCMC) methods. It lets us draw samples from complex probability distributions — even ones we can only evaluate up to an unknown constant — which is exactly what Bayesian inference and many geophysical problems demand. In this post we break the algorithm down and implement it in Python.

What is the Metropolis-Hastings algorithm?

MH produces a sequence of random samples from a target distribution for which direct sampling is hard. It’s especially valuable in high-dimensional spaces or with complex distributions.

Key concepts

1. Target distribution $\pi(x)$: the distribution we want to sample from.
2. Proposal distribution $q(x’ \mid x)$: a simpler distribution used to propose new states.
3. Markov chain: a sequence of states where each depends only on the previous one.

The algorithm

1. Initialization — start with an initial state $x_0$.
2. Proposal step — from the current state $x_t$, propose a new state $x’$ using $q(x’ \mid x_t)$.

3. Acceptance probability — compute

   \[\alpha = \min\left(1, \frac{\pi(x')\, q(x_t \mid x')}{\pi(x_t)\, q(x' \mid x_t)}\right)\]

   This biases the chain toward states with higher probability under the target.

4. Accept/reject — draw $u \sim U(0,1)$. If $u \leq \alpha$, accept ($x_{t+1} = x’$); otherwise reject and stay ($x_{t+1} = x_t$).
5. Iterate — repeat steps 2–4 many times to collect samples.

Python implementation

Example 1: sampling from a Gaussian distribution

import numpy as np
import matplotlib.pyplot as plt

def target_distribution(x):
    return np.exp(-0.5 * x**2)  # Gaussian distribution (unnormalized)

def proposal_distribution(x):
    return np.random.normal(x, 1)  # Propose a new state from a normal distribution centered at x

def metropolis_hastings(initial_state, target_distribution, proposal_distribution, iterations):
    samples = [initial_state]
    current_state = initial_state

    for _ in range(iterations):
        proposed_state = proposal_distribution(current_state)
        acceptance_probability = min(1, target_distribution(proposed_state) / target_distribution(current_state))

        if np.random.rand() < acceptance_probability:
            current_state = proposed_state
        
        samples.append(current_state)
    
    return np.array(samples)

# Parameters
initial_state = 0
iterations = 10000

# Generate samples
samples = metropolis_hastings(initial_state, target_distribution, proposal_distribution, iterations)

# Plotting the results
plt.figure(figsize=(10, 6))
plt.hist(samples, bins=50, density=True, alpha=0.6, color='g', label='Samples')
x = np.linspace(-5, 5, 1000)
plt.plot(x, target_distribution(x), 'r', lw=2, label='Target Distribution')
plt.title('Metropolis-Hastings Sampling from a Gaussian Distribution')
plt.xlabel('x')
plt.ylabel('Density')
plt.legend()
plt.show()

Notice the acceptance line uses only target(proposed) / target(current) — no q terms. That’s because the proposal here is a Gaussian centered on the current state, which is symmetric: $q(x’ \mid x_t) = q(x_t \mid x’)$, so the $q$ factors cancel and the general MH ratio collapses to the simpler Metropolis ratio. It also uses the unnormalized Gaussian — and still recovers the right distribution, exactly as the ratio argument promised.

Example 2: sampling from a multimodal distribution

def target_distribution(x):
    return 0.3 * np.exp(-0.2 * (x - 3)**2) + 0.7 * np.exp(-0.2 * (x + 3)**2)

def proposal_distribution(x):
    return np.random.normal(x, 1)

# Parameters
initial_state = 0
iterations = 10000

# Generate samples
samples = metropolis_hastings(initial_state, target_distribution, proposal_distribution, iterations)

# Plotting the results
plt.figure(figsize=(10, 6))
plt.hist(samples, bins=50, density=True, alpha=0.6, color='g', label='Samples')
x = np.linspace(-10, 10, 1000)
plt.plot(x, target_distribution(x), 'r', lw=2, label='Target Distribution')
plt.title('Metropolis-Hastings Sampling from a Multimodal Distribution')
plt.xlabel('x')
plt.ylabel('Density')
plt.legend()
plt.show()

Recap

Without scrolling up — can you state the loop? Metropolis-Hastings:

• Proposes a new state $x’$ from $q(\cdot \mid x_t)$,
• accepts it with probability $\alpha = \min!\left(1, \frac{\pi(x’)\,q(x_t\mid x’)}{\pi(x_t)\,q(x’\mid x_t)}\right)$, else stays put,
• repeats to build a Markov chain whose stationary distribution is the target $\pi$,
• and works on unnormalized densities because the constant cancels in the ratio.

Experiment with different targets and proposal step sizes to feel how the chain explores — that intuition is what makes the difference in real MCMC applications.

References

1. Equation of State Calculations by Fast Computing Machines — Metropolis, Rosenbluth, Rosenbluth, Teller & Teller, 1953, The Journal of Chemical Physics, 21(6), 1087.
2. Monte Carlo Sampling Methods Using Markov Chains and Their Applications — Hastings, 1970, Biometrika, 57(1), 97–109.

Where to go next

Related posts here: Monte Carlo methods and earthquake location problem and Hypothesis test for the significance of linear trend using Monte Carlo simulations.