The Metropolis–Hastings (MH) algorithm is a Markov Chain Monte Carlo
(MCMC) technique widely used to sample from complex probability
distributions. Specifically, this report employs a Random Walk
Metropolis–Hastings approach, in which proposed updates to the parameter
values are generated by random perturbations around the current state.
This algorithm enables efficient approximation the target distribution
by drawing a series of dependent samples.
Part 1a)
Step 1: Initialise N = 10000, s = 1, x and define the probability
density function
# Initialize total iterations, random walk step size, and initial chain value
total_iter <- 10000
proposal_sd <- 1
initial_val <- 1
# Define the target PDF (Laplace distribution)
target_pdf <- function(u) {
0.5 * exp(-abs(u))}
Step 2: Compute acceptance ratio, sample mean and standard
deviation
# For reproducible results
set.seed(2025)
# Create a vector to store the Markov chain, including the initial state
chain_vec <- numeric(total_iter + 1)
chain_vec[1] <- initial_val
# Generate the chain using a random walk Metropolis–Hastings update
for (idx in 2:(total_iter + 1)) {
# Propose from a Normal centered at the current value
candidate_val <- rnorm(1, mean = chain_vec[idx - 1], sd = proposal_sd)
# Acceptance ratio
ratio_val <- target_pdf(candidate_val) / target_pdf(chain_vec[idx - 1])
if (runif(1) < ratio_val) {
chain_vec[idx] <- candidate_val
} else {
chain_vec[idx] <- chain_vec[idx - 1]
}
}
# Compute sample mean and sample standard deviation
chain_mean <- mean(chain_vec)
chain_sd <- sd(chain_vec)
print(paste("Sample Mean:", chain_mean))
## [1] "Sample Mean: 0.10283493141666"
print(paste("Sample Standard Deviation:", chain_sd))
## [1] "Sample Standard Deviation: 1.37523764167713"
Visualizing Samples vs. True PDF
library(ggplot2)
# Convert the chain into a data frame for plotting
chain_df <- data.frame(samples = chain_vec)
# Histogram of sampled values, overlay with kernel density and the theoretical Laplace PDF
chain_plot <- ggplot(chain_df, aes(x = samples)) +
geom_histogram(aes(y = ..density..), bins = 30, fill = "black", color = "white", alpha = 0.9) +
geom_density(aes(color = "Kernel Density"), show.legend = TRUE) +
stat_function(fun = target_pdf, aes(color = "Laplace PDF"), show.legend = TRUE) +
labs(
title = "Random Walk MH Samples vs. Target PDF",
x = "Sample Values",
y = "Density"
) +
theme_minimal() +
theme(
plot.title = element_text(hjust = 0.5),
plot.background = element_rect(fill = "white")
) +
coord_cartesian(ylim = c(0, 0.6))
print(chain_plot)
## Warning: The dot-dot notation (`..density..`) was deprecated in ggplot2 3.4.0.
## ℹ Please use `after_stat(density)` instead.
## This warning is displayed once every 8 hours.
## Call `lifecycle::last_lifecycle_warnings()` to see where this warning was
## generated.
Part 1b)
Calculate the R-hat for the random walk Metropolis algorithm with N,
s and J fixed
# Define the number of parallel chains(J), chain length(N), and tiny step (S)
num_chains <- 4
chain_length <- 2000
tiny_step <- 0.001
set.seed(2025)
# Initialize vectors to store mean and variance from each chain
multi_means <- numeric(num_chains)
multi_vars <- numeric(num_chains)
for (ch in seq_len(num_chains)) {
# Set up storage for one chain
chain_rw <- numeric(chain_length + 1)
chain_rw[1] <- runif(1, min = 0, max = 10)
# Perform the Metropolis updates
for (i in 2:(chain_length + 1)) {
candidate_point <- rnorm(1, mean = chain_rw[i - 1], sd = tiny_step)
ratio_rw <- target_pdf(candidate_point) / target_pdf(chain_rw[i - 1])
if (runif(1) < ratio_rw) {
chain_rw[i] <- candidate_point
} else {
chain_rw[i] <- chain_rw[i - 1]
}
}
# Compute mean and variance for this chain
m_val <- mean(chain_rw)
v_val <- var(chain_rw)
# Regenerate chain if degenerate
while (m_val == 0 || v_val == 0) {
chain_rw <- numeric(chain_length + 1)
chain_rw[1] <- runif(1, min = 0, max = 10)
for (i in 2:(chain_length + 1)) {
candidate_point <- rnorm(1, mean = chain_rw[i - 1], sd = tiny_step)
ratio_rw <- target_pdf(candidate_point) / target_pdf(chain_rw[i - 1])
if (runif(1) < ratio_rw) {
chain_rw[i] <- candidate_point
} else {
chain_rw[i] <- chain_rw[i - 1]
}
}
m_val <- mean(chain_rw)
v_val <- var(chain_rw)
}
multi_means[ch] <- m_val
multi_vars[ch] <- v_val
}
# Print chain-wise results
for (ch_idx in seq_len(num_chains)) {
print(paste0("Chain ", ch_idx, " Mean: ", multi_means[ch_idx]))
}
## [1] "Chain 1 Mean: 7.31603240758345"
## [1] "Chain 2 Mean: 1.36686917489894"
## [1] "Chain 3 Mean: 8.39301887503724"
## [1] "Chain 4 Mean: 0.62088159832511"
for (ch_idx in seq_len(num_chains)) {
print(paste0("Chain ", ch_idx, " Variance: ", multi_vars[ch_idx]))
}
## [1] "Chain 1 Variance: 0.000241227451640372"
## [1] "Chain 2 Variance: 0.000379619349923515"
## [1] "Chain 3 Variance: 0.000150828610372669"
## [1] "Chain 4 Variance: 0.000356809962395742"
Checking if algorithm has converged using R-hat value
within_chain_var <- mean(multi_vars)
grand_mean <- mean(multi_means)
between_chain_var <- sum((multi_means - grand_mean)^2) / num_chains
r_hat_val <- sqrt((between_chain_var + within_chain_var) / within_chain_var)
print(paste("R-hat value:", r_hat_val))
## [1] "R-hat value: 206.084819413271"
print(paste("Within-chain variance (W):", within_chain_var))
## [1] "Within-chain variance (W): 0.000282121343583074"
print(paste("Grand mean (M):", grand_mean))
## [1] "Grand mean (M): 4.42420051396118"
print(paste("Between-chain variance (B):", between_chain_var))
## [1] "Between-chain variance (B): 11.9816801437582"
Keeping N and J fixed, repeat the algorithm and compute R-hat values
with different s values
set.seed(2025)
num_experiments <- 50
store_s_vals <- numeric(num_experiments)
store_r_hat <- numeric(num_experiments)
for (exp_idx in 1:num_experiments) {
# Draw a random step size from [0.001, 1]
rand_s <- runif(1, min = 0.001, max = 1)
store_s_vals[exp_idx] <- rand_s
# Prepare arrays for chain stats
chain_means_temp <- numeric(num_chains)
chain_vars_temp <- numeric(num_chains)
# Generate all chains for this step size
for (j in seq_len(num_chains)) {
local_chain <- numeric(chain_length + 1)
local_chain[1] <- runif(1, min = 1, max = 10)
for (i in 2:(chain_length + 1)) {
prop_step <- rnorm(1, mean = local_chain[i - 1], sd = rand_s)
ratio_step <- target_pdf(prop_step) / target_pdf(local_chain[i - 1])
if (runif(1) < ratio_step) {
local_chain[i] <- prop_step
} else {
local_chain[i] <- local_chain[i - 1]
}
}
temp_m <- mean(local_chain)
temp_v <- var(local_chain)
# Re-generate if degenerate
while (temp_m == 0 || temp_v == 0) {
local_chain <- numeric(chain_length + 1)
local_chain[1] <- runif(1, min = 1, max = 10)
for (i in 2:(chain_length + 1)) {
prop_step <- rnorm(1, mean = local_chain[i - 1], sd = rand_s)
ratio_step <- target_pdf(prop_step) / target_pdf(local_chain[i - 1])
if (runif(1) < ratio_step) {
local_chain[i] <- prop_step
} else {
local_chain[i] <- local_chain[i - 1]
}
}
temp_m <- mean(local_chain)
temp_v <- var(local_chain)
}
chain_means_temp[j] <- temp_m
chain_vars_temp[j] <- temp_v
}
# Compute R-hat for the current step size
w_val <- mean(chain_vars_temp)
m_val <- mean(chain_means_temp)
b_val <- sum((chain_means_temp - m_val)^2) / num_chains
store_r_hat[exp_idx] <- sqrt((b_val + w_val) / w_val)
}
# Print out the R-hat results
for (i in seq_along(store_r_hat)) {
cat("R-hat value for experiment", i, ":", store_r_hat[i], "\n")
}
## R-hat value for experiment 1 : 1.001214
## R-hat value for experiment 2 : 1.003216
## R-hat value for experiment 3 : 1.003943
## R-hat value for experiment 4 : 1.006484
## R-hat value for experiment 5 : 1.006643
## R-hat value for experiment 6 : 1.014261
## R-hat value for experiment 7 : 1.003671
## R-hat value for experiment 8 : 1.001425
## R-hat value for experiment 9 : 1.003522
## R-hat value for experiment 10 : 1.008971
## R-hat value for experiment 11 : 1.008042
## R-hat value for experiment 12 : 1.025151
## R-hat value for experiment 13 : 1.010688
## R-hat value for experiment 14 : 1.058981
## R-hat value for experiment 15 : 1.005817
## R-hat value for experiment 16 : 1.018423
## R-hat value for experiment 17 : 1.011306
## R-hat value for experiment 18 : 1.019429
## R-hat value for experiment 19 : 1.063338
## R-hat value for experiment 20 : 1.015571
## R-hat value for experiment 21 : 1.006284
## R-hat value for experiment 22 : 1.033345
## R-hat value for experiment 23 : 1.003507
## R-hat value for experiment 24 : 1.002972
## R-hat value for experiment 25 : 1.009715
## R-hat value for experiment 26 : 1.005645
## R-hat value for experiment 27 : 1.002291
## R-hat value for experiment 28 : 1.001018
## R-hat value for experiment 29 : 1.003177
## R-hat value for experiment 30 : 1.003575
## R-hat value for experiment 31 : 1.11584
## R-hat value for experiment 32 : 1.011023
## R-hat value for experiment 33 : 1.002769
## R-hat value for experiment 34 : 1.017684
## R-hat value for experiment 35 : 1.001204
## R-hat value for experiment 36 : 1.329461
## R-hat value for experiment 37 : 1.009586
## R-hat value for experiment 38 : 1.025362
## R-hat value for experiment 39 : 1.063411
## R-hat value for experiment 40 : 1.041186
## R-hat value for experiment 41 : 1.046596
## R-hat value for experiment 42 : 1.004053
## R-hat value for experiment 43 : 1.026688
## R-hat value for experiment 44 : 1.00818
## R-hat value for experiment 45 : 1.002364
## R-hat value for experiment 46 : 1.514574
## R-hat value for experiment 47 : 1.090423
## R-hat value for experiment 48 : 1.001141
## R-hat value for experiment 49 : 1.003697
## R-hat value for experiment 50 : 1.005901
Plot r-hat values over s values
# Combine s and R-hat into a data frame for ggplot
results_df <- data.frame(
step_sizes = store_s_vals,
r_hat_vals = store_r_hat
)
rhat_plot <- ggplot(results_df, aes(x = step_sizes, y = r_hat_vals)) +
geom_point() +
geom_smooth(method = "lm", se = FALSE, color = "red") +
labs(
title = "R-hat vs. Step Size",
x = "Step Size (s)",
y = "R-hat Values"
) +
theme_minimal() +
theme(
plot.title = element_text(hjust = 0.5),
plot.background = element_rect(fill = "white")
) +
coord_cartesian(ylim = c(0, 2)) +
geom_hline(yintercept = 1.05, linetype = "dashed", color = "blue")
print(rhat_plot)
## `geom_smooth()` using formula = 'y ~ x'
