Metropolis-Hastings library(ggplot2) library(tidyverse) library(gridExtra) library(MASS) Random Walk Metropolis Hastings Let’s implement a Random Walk Metropolis-Hastings algorithm with a symmetric proposal density.
rwmh <- function(start, niter, proposal, target){ # Store accepted, rejected and all samples accepted <- matrix(NA, nrow=niter, dimnames=list(NULL, "accepted")) rejected <- matrix(NA, nrow=niter, dimnames=list(NULL, "rejected")) samples <- matrix(NA, nrow=niter) # Set first sample to start and sample uniform random numbers for efficiency z <- start accepted[1] <- z; rejected[1] <- z; samples[1] <- z u <- runif(niter) # Draw candidate, calculate acceptance prob for niter-1 times for (i in 2:niter){ cand <- proposal$sample(z) prob <- min(1, target$density(cand)/target$density(z)) if (u[i] <= prob) {accepted[i] <- cand; z <- cand} else {rejected[i] <- cand} samples[i] <- z } result <- list("samples"=samples, "accepted"=accepted, "rejected"=rejected, "proposal"=proposal, "target"=target) return(result) } Plotting Function The following function takes the output of the Random-Walk Metropolis-Hastings algorithm and it produces two plots:
