Metropolis-Hastings in R

library(purrr)
library(ggplot2)
library(stringr)

Specifying the Target and Proposal

# density of the target distribution
dtarget <- partial(dchisq, df=10) 

# draws a proposal from uniform +-1 around the previous state
propose <- function(prev){runif(1, min=prev-1, max=prev+1)}  

# the density of the proposal (conditional on the previous state)
dproposal <- function(prev, proposal){1/2}

# a boolean that represents rejection/acceptance of the proposal
accept <- function(prev, proposed){
  U <- runif(1)
  accept_prob <- min(1, dtarget(proposed) / dtarget(prev) * 
                        dproposal(proposed, prev) / dproposal(prev, proposed))
  return(U < accept_prob)
}

Running Metropolis-Hastings

# initializing
n <- 100000  # change it if your computer can't handle 1e5 samples.
previous_state <- 1
samples <- numeric(length = n)

# M-H
for (k in 1:n){
  proposal <- propose(previous_state)
  accepted <- accept(previous_state, proposal)
  if (accepted){
    previous_state <- proposal
    samples[k] <- proposal
  } else {
    samples[k] <- previous_state
  }
}

Plotting the Samples

binwidth <- 0.5  # how fine the density/bar plot should be

ggplot(data = data.frame(samples=samples)) +
  geom_histogram(aes(x=samples,
                    y = after_stat(count / sum(count) / binwidth),
                    fill="M-H Samples"),
                 binwidth = binwidth,
                 alpha=0.6) +
  geom_function(aes(color="Target Distribution"),
                fun = dtarget) +
  labs(x="Sample Value", y="Relative Density") +
  scale_fill_manual(name = "", values = c("M-H Samples" = "#233F7E")) +
  scale_color_manual(name = "", values = c("Target Distribution" = "#233F7E")) +
  ggtitle(str_interp("Distribution of ${n} M-H Samples")) +
  theme_minimal() + theme(legend.position="bottom")