MH Algo

library(ggplot2)
library(patchwork)  
library(mvtnorm)  
set.seed(431)

One Metropolis–Hastings update step

MHstep <- function(theta_current, sig) {
  theta_star <- rnorm(1, mean = theta_current, sd = sig)  # propose candidate θ*
  accprob <- targetdist(theta_star) / targetdist(theta_current)  # acceptance probability
  accprob <- min(1, accprob)
  # accept or reject
  u <- runif(1)
  if (u <= accprob) {
    theta_next <- theta_star
    a <- 1
  } else {
    theta_next <- theta_current
    a <- 0
  }
  return(list(theta1 = theta_next, a = a))
}

Metropolis–Hastings Sampler

MH_sampler <- function(nsamp, sig, burnin, lag, theta_current) {
  theta_samples <- numeric(nsamp)
  acc <- c(accepted = 0, proposed = 0)
  theta <- theta_current   # starting value
  # Burn-in
  for (i in 1:burnin) {
    step <- MHstep(theta, sig)
    theta <- step$theta1
    acc <- acc + c(step$a, 1)
  }
  # Main sampling
  for (i in 1:nsamp) {
    for (j in 1:lag) {
      step <- MHstep(theta, sig)
      theta <- step$theta1
      acc <- acc + c(step$a, 1)
    }
    theta_samples[i] <- theta
  }
  return(list(samples = theta_samples, acc = acc))
}

Visualization function

plot_MH <- function(sig, nsamp = NULL, theta_initial = NULL) {
  
  # Pull global defaults if user does not override
  if (is.null(nsamp)) nsamp <- get("nsamp", envir = .GlobalEnv)
  if (is.null(theta_initial)) theta_initial <- get("theta_initial", envir = .GlobalEnv)
  
  # Run MH sampler
  result <- MH_sampler(
    nsamp = nsamp,
    sig = sig,
    burnin = burnin,
    lag = lag,
    theta_current = theta_initial   # pass renamed argument to MH sampler
  )
  
  theta_vals <- result$samples
  
  # True normalized density
  thetas <- seq(-3, 3, length.out = 500)
  dens <- targetdist(thetas)
  dens <- dens / sum(dens) / (thetas[2] - thetas[1])
  
  df_samples <- data.frame(theta = theta_vals)
  df_true    <- data.frame(theta = thetas, y = dens)
  
  # Histogram + True Density 
 p1 <- ggplot() +
  geom_histogram(data = df_samples,
                 aes(theta, y = ..density.., fill = "Samples"),
                 bins = 30, color = "black") +
  geom_line(data = df_true,
            aes(theta, y, color = "True Density"),
            size = 1.1) +
  scale_color_manual(values = c("True Density" = "red")) +
  scale_fill_manual(values = c("Samples" = "gray80")) +
  labs(
    title = paste0("Metropolis–Hastings Sampling (σ = ", sig, ")"),
    subtitle = paste0(
      "burn-in = ", burnin,
      ",  nsamp = ", nsamp,
      ",  initial θ = ", theta_initial
    ),
    x = expression(theta), y = "density",
    color = "", fill = ""
  ) +
  theme_bw(base_size = 14)

  
  # ---------------- Trace Plot ----------------
  df_trace <- data.frame(iter = 1:length(theta_vals), value = theta_vals)
  
  p2 <- ggplot(df_trace, aes(iter, value)) +
    geom_line(color = "black") +
    labs(
      title = "Trace Plot",
      subtitle = paste0(
        "burn-in = ", burnin,
        ",  kept samples = ", nsamp,
        ",  initial θ = ", theta_initial
      ),
      y = expression(theta~value),
      x = "iteration"
    ) +
    theme_bw(base_size = 14)
  
  # Return stacked plots
  p1 / p2
}

Example

set.seed(431)
#   Target distribution (unnormalized)
targetdist <- function(theta) {
  exp(-theta^4/6 + theta^2/2) * cos(theta^3/4)
}
burnin <- 8000
lag <- 1
nsamp <- 50000
theta_initial <- 0   

# global default
#plot_MH(sig = 0.1)
#plot_MH(sig = 1)
#plot_MH(sig = 10)

# 1. Well-tuned, centered start
plot_MH(sig = 1, nsamp = 3000, theta_initial = 0)

# 2. Same σ, bad initial value 
plot_MH(sig = 1, nsamp = 3000, theta_initial = 6)

# 3. Smaller step size, slower mixing
plot_MH(sig = 0.01, nsamp = 3000, theta_initial = 0)

Gibbs Sampler for Normal Model with Semi-Conjugate Priors

# Data summary
n    <- 24
ybar <- 7.8730           # sample mean
s     <- 0.05353         # sample sd

### Prior parameters
mu0       <- 7.9876      # prior mean for mu
sigma0.2  <- 0.000625    # prior variance of mu
alpha     <- 2.5         # IG(alpha, beta)
beta      <- 0.0043

#  Run Gibbs sampler
Nsim <- 100000               # total draws
mus <- numeric(Nsim)         # store mu draws
sigma2 <- numeric(Nsim)      # store sigma^2 draws

## Initial values
mus[1] <- 8
SSE1   <- (n-1)*s^2 + n*(ybar - mus[1])^2
sigma2[1] <- 1 / rgamma(1, n/2 + alpha, SSE1/2 + beta)

## Gibbs updates
for (k in 2:Nsim) {
  
  ### 1. Sample mu | sigma^2, y
  var_mu  <- 1 / (n / sigma2[k-1] + 1 / sigma0.2)
  mean_mu <- var_mu * (n*ybar / sigma2[k-1] + mu0 / sigma0.2)
  mus[k] <- rnorm(1, mean_mu, sqrt(var_mu))
  
  ### 2. Sample sigma^2 | mu, y
  SSE <- (n-1)*s^2 + n*(ybar - mus[k])^2
  sigma2[k] <- 1 / rgamma(1, n/2 + alpha, SSE/2 + beta)
}


#  Posterior summaries
cat("Posterior mean(mu):", mean(mus), "\n")

Posterior mean(mu): 7.893102 

cat("Posterior sd(mu):",   sd(mus), "\n")

Posterior sd(mu): 0.01170627 

cat("95% CI for mu:\n"); print(quantile(mus, c(0.025, 0.975)))

95% CI for mu:
    2.5%    97.5% 
7.871854 7.917967 

cat("\nPosterior mean(sigma^2):", mean(sigma2), "\n")

Posterior mean(sigma^2): 0.003244212 

cat("95% CI for sigma^2:\n"); print(quantile(sigma2, c(0.025, 0.975)))

95% CI for sigma^2:
       2.5%       97.5% 
0.001803182 0.005834502 

#  Visualize first few Gibbs steps
maxit <- 20
plot(mus[1:maxit], sigma2[1:maxit], type = "n",
     xlab = expression(mu), ylab = expression(sigma^2),
     main = "First 20 Gibbs Iterations")
arrows(mus[1:(maxit-1)], sigma2[1:(maxit-1)],
       mus[2:maxit], sigma2[2:maxit],
       length = 0.1, angle = 30, col = rainbow(maxit-1))

points(mus[1:maxit], sigma2[1:maxit], pch = 20)

#  Trace plots
trace_mu <- ggplot(data.frame(iter = 1:Nsim, mu = mus),
                   aes(iter, mu)) +
  geom_line(color = "steelblue") +
  labs(title = "Trace Plot for μ",
       x = "Iteration", y = expression(mu)) +
  theme_bw(14)
trace_sigma2 <- ggplot(data.frame(iter = 1:Nsim, sigma2 = sigma2),
                       aes(iter, sigma2)) +
  geom_line(color = "firebrick") +
  labs(title = expression("Trace Plot for " * sigma^2),
       x = "Iteration", y = expression(sigma^2)) +
  theme_bw(14)

#  Posterior density plots
dens_mu <- ggplot(data.frame(mu = mus), aes(mu)) +
  geom_density(fill = "steelblue", alpha = 0.4) +
  geom_vline(xintercept = mean(mus), color = "blue", linetype = 2) +
  labs(title = "Posterior Density of μ",
       x = expression(mu), y = "Density") +
  theme_bw(14)
dens_sigma2 <- ggplot(data.frame(sigma2 = sigma2), aes(sigma2)) +
  geom_density(fill = "firebrick", alpha = 0.4) +
  geom_vline(xintercept = mean(sigma2), color = "red", linetype = 2) +
  labs(title = "Posterior Density of σ²",
       x = expression(sigma^2), y = "Density") +
  theme_bw(14)

#  Joint posterior (μ, σ²)
df_post <- data.frame(mu = mus, sigma2 = sigma2)
joint_contour <- ggplot(df_post, aes(mu, sigma2)) +
  stat_density_2d(aes(fill = ..level..), geom = "polygon", alpha = 0.6) +
  scale_fill_viridis_c() +
  labs(title = "Posterior Joint Density of (μ, σ²)",
       x = expression(mu), y = expression(sigma^2)) +
  theme_bw(14)
joint_heatmap <- ggplot(df_post, aes(mu, sigma2)) +
  stat_density_2d_filled(contour_var = "density") +
  labs(title = "Posterior Joint Density Heatmap",
       x = expression(mu), y = expression(sigma^2)) +
  theme_bw(14)


#  Plot combined panels
( trace_mu | trace_sigma2 ) / ( dens_mu | dens_sigma2 )

joint_contour

joint_heatmap

MH Sampler

set.seed(431)
## Log-posterior for (mu, log_sigma2)
log_post <- function(mu, log_sigma2) {
  sigma2 <- exp(log_sigma2)
  
  # SSE(mu) using sufficient stats
  SSE <- (n - 1)*s^2 + n*(ybar - mu)^2
  
  # log-likelihood: N(mu, sigma2) for n observations
  lp_y <- -(n/2) * log(2*pi*sigma2) - SSE/(2*sigma2)
  
  # prior for mu ~ N(mu0, sigma0^2)
  lp_mu <- dnorm(mu, mean = mu0, sd = sqrt(sigma0.2), log = TRUE)
  
  # prior for sigma2 ~ IG(alpha, beta), up to constants:
  # p(sigma2) ∝ (sigma2)^-(alpha+1) exp(-beta/sigma2)
  lp_s2 <- -(alpha + 1)*log(sigma2) - beta/sigma2
  
  # Jacobian for transform sigma2 = exp(log_sigma2)
  lp_jac <- log_sigma2
  return(lp_y + lp_mu + lp_s2 + lp_jac)
} 
  
  
mh_sampler <- function(Nsim, start_mu, start_log_s2,
                       prop_sd_mu, prop_sd_log_s2) {
  mu     <- numeric(Nsim)
  log_s2 <- numeric(Nsim)
  mu[1]     <- start_mu
  log_s2[1] <- start_log_s2
  acc <- 0
  
  for (k in 2:Nsim) {
    # Propose new (mu, log_sigma2) with Gaussian random walk
    mu_prop     <- rnorm(1, mean = mu[k-1],     sd = prop_sd_mu)
    log_s2_prop <- rnorm(1, mean = log_s2[k-1], sd = prop_sd_log_s2)
    
    # Log acceptance ratio (symmetric proposal ⇒ q cancels)
    log_alpha <- log_post(mu_prop, log_s2_prop) -
                 log_post(mu[k-1], log_s2[k-1])
    
    if (log(runif(1)) < log_alpha) {
      # accept
      mu[k]     <- mu_prop
      log_s2[k] <- log_s2_prop
      acc <- acc + 1
    } else {
      # reject (stay at previous)
      mu[k]     <- mu[k-1]
      log_s2[k] <- log_s2[k-1]
    }
  }
  
  list(mu = mu,
       sigma2 = exp(log_s2),
       acc_rate = acc / (Nsim - 1))
}

# Use Gibbs initial sigma2[1] as starting point for MH
mh_out <- mh_sampler(
  Nsim          = Nsim,
  start_mu      = mus[1],
  start_log_s2  = log(sigma2[1]),
  prop_sd_mu    = 0.005,   # tune these
  prop_sd_log_s2 = 0.05
)

mh_out$acc_rate   # check acceptance rate; aim ~ 0.2–0.5

[1] 0.8239982

Comparsion Plots

library(ggplot2)
library(patchwork)

df_comp <- data.frame(
  iter          = 1:Nsim,
  gibbs_mu      = mus,
  mh_mu         = mh_out$mu,
  gibbs_sigma2  = sigma2,
  mh_sigma2     = mh_out$sigma2
)

# Trace Plot
trace_mu_both <- ggplot(df_comp) +
  geom_line(aes(iter, gibbs_mu, colour = "Gibbs")) +
  geom_line(aes(iter, mh_mu,    colour = "MH"), alpha = 0.7) +
  scale_colour_manual(values = c("Gibbs" = "steelblue",
                                 "MH"    = "orange")) +
  labs(title = "Trace Plot for μ: Gibbs vs MH",
       x = "Iteration", y = expression(mu), colour = "") +
  theme_bw(14)
trace_s2_both <- ggplot(df_comp) +
  geom_line(aes(iter, gibbs_sigma2, colour = "Gibbs")) +
  geom_line(aes(iter, mh_sigma2,    colour = "MH"), alpha = 0.7) +
  scale_colour_manual(values = c("Gibbs" = "firebrick",
                                 "MH"    = "darkgreen")) +
  labs(title = expression("Trace Plot for " * sigma^2 * ": Gibbs vs MH"),
       x = "Iteration", y = expression(sigma^2), colour = "") +
  theme_bw(14)
trace_mu_both / trace_s2_both

# Density comparison
dens_mu_both <- ggplot(df_comp) +
  geom_density(aes(gibbs_mu, colour = "Gibbs"), adjust = 1.2) +
  geom_density(aes(mh_mu,    colour = "MH"),    adjust = 1.2, linetype = 2) +
  scale_colour_manual(values = c("Gibbs" = "steelblue",
                                 "MH"    = "orange")) +
  labs(title = "Posterior of μ: Gibbs vs MH",
       x = expression(mu), y = "Density", colour = "") +
  theme_bw(14)
dens_s2_both <- ggplot(df_comp) +
  geom_density(aes(gibbs_sigma2, colour = "Gibbs"), adjust = 1.2) +
  geom_density(aes(mh_sigma2,    colour = "MH"),    adjust = 1.2, linetype = 2) +
  scale_colour_manual(values = c("Gibbs" = "firebrick",
                                 "MH"    = "darkgreen")) +
  labs(title = expression("Posterior of " * sigma^2 * ": Gibbs vs MH"),
       x = expression(sigma^2), y = "Density", colour = "") +
  theme_bw(14)
dens_mu_both | dens_s2_both

# --- Gibbs heatmap 
mh_mu_q  <- quantile(df_comp$mh_mu,     probs = c(0.005, 0.995))
mh_s2_q  <- quantile(df_comp$mh_sigma2, probs = c(0.005, 0.995))
zoom_x <- c(mh_mu_q[1], mh_mu_q[2])
zoom_y <- c(mh_s2_q[1], mh_s2_q[2])
heat_gibbs_zoom <- ggplot(df_comp, aes(gibbs_mu, gibbs_sigma2)) +
  stat_density_2d_filled(contour_var = "density", alpha = 0.9) +
  stat_density_2d(color = "white", linewidth = 0.4) +
  scale_fill_viridis_d(option = "C") +
  coord_cartesian(xlim = zoom_x, ylim = zoom_y, expand = FALSE) +
  labs(
    title = "Joint Posterior Heatmap (Gibbs)",
    x = expression(mu),
    y = expression(sigma^2),
    fill = "Density Level"
  ) +
  theme_bw(base_size = 16) +
  theme(
    plot.title = element_text(size = 20, face = "bold", hjust = 0.5),
    axis.title = element_text(size = 18),
    axis.text  = element_text(size = 14),
    legend.title = element_text(size = 14),
    legend.text  = element_text(size = 12),
    panel.border = element_rect(linewidth = 0.8, color = "black"),
    panel.grid = element_blank()
  )
heat_gibbs_zoom

# MH heatmap 

heat_mh_zoom <- ggplot(df_comp, aes(mh_mu, mh_sigma2)) +
  stat_density_2d_filled(contour_var = "density", alpha = 0.9) +
  stat_density_2d(color = "white", linewidth = 0.4) +
  scale_fill_viridis_d(option = "C") +
  coord_cartesian(xlim = zoom_x, ylim = zoom_y, expand = FALSE) +
  labs(
    title = "Joint Posterior Heatmap (MH)",
    x = expression(mu),
    y = expression(sigma^2),
    fill = "Density Level"
  ) +
  theme_bw(base_size = 16) +
  theme(
    plot.title = element_text(size = 20, face = "bold", hjust = 0.5),
    axis.title = element_text(size = 18),
    axis.text  = element_text(size = 14),
    legend.title = element_text(size = 14),
    legend.text  = element_text(size = 12),
    panel.border = element_rect(linewidth = 0.8, color = "black"),
    panel.grid = element_blank()
  )
heat_mh_zoom

# Set layout: 1 row, 2 columns
par(mfrow = c(1, 2), mar = c(5, 5, 4, 2))

# ------------------------------
# Panel 1: Posterior of μ
# ------------------------------
dens_gibbs_mu <- density(df_comp$gibbs_mu)
dens_mh_mu    <- density(df_comp$mh_mu)

plot(dens_gibbs_mu,
     main = "Posterior of μ: Gibbs vs MH",
     xlab = expression(mu),
     ylab = "Density",
     col  = "steelblue",
     lwd  = 2,
     lty  = 1,          # solid
     cex.lab = 1.3,
     cex.main = 1.4)

lines(dens_mh_mu,
      col = "orange",
      lwd = 2,
      lty = 2)          # dashed

legend("topright",
       legend = c("Gibbs (solid)", "MH (dashed)"),
       col = c("steelblue", "orange"),
       lwd = 2,
       lty = c(1, 2),
       bty = "n",
       cex = 1.1)

# ------------------------------
# Panel 2: Posterior of σ²
# ------------------------------
dens_gibbs_s2 <- density(df_comp$gibbs_sigma2)
dens_mh_s2    <- density(df_comp$mh_sigma2)

plot(dens_gibbs_s2,
     main = expression("Posterior of " * sigma^2 * ": Gibbs vs MH"),
     xlab = expression(sigma^2),
     ylab = "Density",
     col  = "firebrick",
     lwd  = 2,
     lty  = 1,          # solid
     cex.lab = 1.3,
     cex.main = 1.4)

lines(dens_mh_s2,
      col = "darkgreen",
      lwd = 2,
      lty = 2)          # dashed

legend("topright",
       legend = c("Gibbs (solid)", "MH (dashed)"),
       col = c("firebrick", "darkgreen"),
       lwd = 2,
       lty = c(1, 2),
       bty = "n",
       cex = 1.1)

# Restore default layout
par(mfrow = c(1, 1))