Random Walk, Metropolis-Hastings Sampling

Introduction (Qs. 1)

The Metropolis-Hastings algorithm uses Markov Chain Monte Carlo method to sample from distributions where direct sampling is difficult.

Algorithm Description

The Metropolis-Hastings algorithm consists of the following steps:

Initialize \(x^{(0)} \sim q(x)\).

For each iteration \(i = 1, 2, \dots\)

• Propose \(x^{\text{cand}} \sim q(x | x^{(i-1)})\).

• Calculate acceptance probability \(\alpha\):

\[ \alpha(x^{\text{cand}} | x^{(i-1)}) = \min \left\{ 1, \frac{q(x^{(i-1)} | x^{\text{cand}}) \pi(x^{\text{cand}})}{q(x^{\text{cand}} | x^{(i-1)}) \pi(x^{(i-1)})} \right\} \]

• Draw \(u \sim \text{Uniform}(0, 1)\).

• If \(u < \alpha\), accept the proposal: \(x^{(i)} \gets x^{\text{cand}}\).

• Otherwise, reject the proposal: \(x^{(i)} \gets x^{(i-1)}\).

Repeat this process for the desired number of iterations.

Implementation in R

# Metropolis-Hastings sampler function
mh_sampler <- function(n, init, target, proposal) {
  samples <- numeric(n)
  current <- init
  current_prob <- target(current)
  
  for (i in 1:n) {
    proposed <- proposal(current)
    proposed_prob <- target(proposed)
    
    acceptance_prob <- proposed_prob / current_prob
    
    if (runif(1) < acceptance_prob) {
      current <- proposed
      current_prob <- proposed_prob
    }
    samples[i] <- current
  }
  samples
}

# Target distribution: N(0, 1)
target_distribution <- function(x) {
  dnorm(x, mean = 0, sd = 1)
}

# Proposal distribution: N(current, 0.5)
proposal_distribution <- function(current) {
  rnorm(1, mean = current, sd = 0.5)
}

# Running the sampler
set.seed(786)   #For reproducibility
sample_Z <- mh_sampler(10000, init = 0, target = target_distribution, proposal = proposal_distribution)

Results

We plot the histogram of the samples obtained via the Metropolis-Hastings algorithm overlaid with the true density function.

# Plotting the results
hist(sample_Z, probability = TRUE, breaks = 40, main = "Samples from MH Sampler from Standard Normal Distribution")
curve(dnorm(x, mean = 0, sd = 1), add = TRUE, col = 'black', lwd = 2)

---

Introduction (Qs. 2)

The RWM algorithm is a special case of the Metropolis-Hastings algorithm where the proposal distribution is centered around the current position and is typically symmetric.

Algorithm Description

The Random Walk Metropolis algorithm consists of the following steps:

For each iteration \(i = 1, 2, \dots\)

• Generate a candidate sample \(x^{\text{cand}}\) from a proposal distribution \(q(x^{(i-1)}, x)\) centered at the current sample \(x^{(i-1)}\). Commonly, \(q\) is a normal distribution, making this a random walk.

• Calculate the acceptance probability \(\alpha\):

\[ \alpha(x^{\text{cand}} | x^{(i-1)}) = \min \left\{ 1, \frac{\pi(x^{\text{cand}})}{\pi(x^{(i-1)})} \right\} \]

• Draw \(u \sim \text{Uniform}(0, 1)\).

• If \(u < \alpha\), accept the proposal: \(x^{(i)} \gets x^{\text{cand}}\).

• Otherwise, reject the proposal: \(x^{(i)} \gets x^{(i-1)}\).

Repeat this process for the desired number of iterations.

R code implementation for Laplace Distribution:

Let’s now translate this algorithm into R code, using the Laplace distribution as the target distribution. The Laplace distribution, with location parameter \(\mu\) and scale parameter \(b\), has the density function: \[ f(x | \mu, b) = \frac{1}{2b} \exp\left(-\frac{|x - \mu|}{b}\right) \] Here’s a basic implementation in R:

# Random Walk Metropolis sampler for Laplace distribution
rwm_laplace <- function(n, init, mu, b, proposal_sd) {
  samples <- numeric(n)
  samples[1] <- init
  target <- function(x) {
    # Laplace distribution density function
    0.5 / b * exp(-abs(x - mu) / b)
  }
  for (i in 2:n) {
    current <- samples[i - 1]
    proposed <- rnorm(1, mean = current, sd = proposal_sd)
    acceptance_prob <- min(1, target(proposed) / target(current))
    
    if (runif(1) < acceptance_prob) {
      samples[i] <- proposed
    } else {
      samples[i] <- current
    }
  }
  samples
}

# Example usage:
set.seed(786)  # For reproducibility
mu <- 0        # Location parameter for Laplace
b <- 1         # Scale parameter for Laplace
n <- 10000
initial_value <- 0
proposal_sd <- 2  # Standard deviation of the proposal distribution (tune this as needed)

sample_L <- rwm_laplace(n, initial_value, mu, b, proposal_sd)

Results

We plot the histogram of the samples obtained via the Random Walk Metropolis algorithm overlaid with the true density function.

# Plotting the results
hist(sample_L, probability = TRUE, breaks = 40, main = "Samples from RWM Sampler for Laplace Distribution")
curve(dexp(abs(x - mu) / b) / b, add = TRUE, col = 'black', lwd = 2)  # Adding theoretical density