Understanding the MetropolisHastings Algorithm: A StepbyStep Guide with Python Examples
The MetropolisHastings algorithm is a cornerstone of Markov Chain Monte Carlo (MCMC) methods, enabling us to generate samples from complex probability distributions. This post breaks down the algorithm and provides Python examples to illustrate its use.
The MetropolisHastings algorithm is a cornerstone of Markov Chain Monte Carlo (MCMC) methods. It allows us to generate samples from complex probability distributions, which is particularly useful in Bayesian statistics and many other fields. In this post, we’ll break down the algorithm and provide Python examples to illustrate its use.
Similar posts
What is the MetropolisHastings Algorithm?
The MetropolisHastings algorithm is used to obtain a sequence of random samples from a target probability distribution for which direct sampling is challenging. It’s particularly valuable when dealing with highdimensional spaces or complex distributions.
Key Concepts
 **Target Distribution $\pi (x)$ **: The probability distribution we want to sample from.
 **Proposal Distribution $q(x’ \vert x)$ **: A simpler distribution used to propose new states in the Markov chain.
 Markov Chain: A sequence of states where the probability of each state depends only on the previous state.
The Algorithm
Here are the steps involved in the MetropolisHastings algorithm:
 Initialization:
 Start with an initial state $x_0$.
 Proposal Step:
 From the current state $ x_t $, propose a new state $ x’ $ using the proposal distribution $q(x’ \vert x_t)$.
 Acceptance Probability:

Calculate the acceptance probability $\alpha$: \begin{equation} \begin{split} \alpha = \min\left(1, \frac{\pi(x’) q(x_t \vert x’)}{\pi(x_t) q(x’ \vert x_t)}\right) \end{split} \end{equation}

This ensures the chain moves towards states with higher probability according to the target distribution.

 Acceptance/Rejection:
 Generate a uniform random number $ u $ from $ U(0, 1) $.
 If $ u \leq \alpha $, accept the proposed state $ x’ $ (i.e., set $ x_{t+1} = x’ $).
 Otherwise, reject $ x’ $ and stay at the current state (i.e., set $ x_{t+1} = x_t $).
 Iteration:
 Repeat steps 2 to 4 for a large number of iterations to obtain a sequence of samples.
Python Implementation
Let’s see how we can implement the MetropolisHastings algorithm in Python.
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('MetropolisHastings Sampling from a Gaussian Distribution')
plt.xlabel('x')
plt.ylabel('Density')
plt.legend()
plt.show()
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('MetropolisHastings Sampling from a Multimodal Distribution')
plt.xlabel('x')
plt.ylabel('Density')
plt.legend()
plt.show()
Conclusion
The MetropolisHastings algorithm is a powerful tool for sampling from complex distributions. By following the steps outlined above and using the provided Python examples, you can implement this algorithm to tackle a variety of problems in statistics and data science. Experiment with different target distributions and proposal mechanisms to see how the algorithm performs in different scenarios.
Disclaimer of liability
The information provided by the Earth Inversion is made available for educational purposes only.
Whilst we endeavor to keep the information uptodate and correct. Earth Inversion makes no representations or warranties of any kind, express or implied about the completeness, accuracy, reliability, suitability or availability with respect to the website or the information, products, services or related graphics content on the website for any purpose.
UNDER NO CIRCUMSTANCE SHALL WE HAVE ANY LIABILITY TO YOU FOR ANY LOSS OR DAMAGE OF ANY KIND INCURRED AS A RESULT OF THE USE OF THE SITE OR RELIANCE ON ANY INFORMATION PROVIDED ON THE SITE. ANY RELIANCE YOU PLACED ON SUCH MATERIAL IS THEREFORE STRICTLY AT YOUR OWN RISK.
Leave a comment