Metropolis-Hastings algorithm in R: correct results?

markov-chains, r

Solution

You had 3 issues:

1) You seem to want to simulate a binomial distribution, so your random walk should be over integers in the range `1:n` rather than real numbers in the range `[0,1]`.

2) You had the numerator and the denominator switched in the computation of `alpha`

3) You had a typo in `X[t] < X[t - 1]`.

Fixing these and cleaning up your code a bit (including using the `dbinom` function as @ZheyuanLi suggested) yields:

metropolisAlgorithm <- function(n, p, T){
  # implementation of the algorithm
  # @n,p binomial parameters
  # @T number of time steps
  X <- rep(0,T)
  X[1] <- sample(1:n,1)
  for (t in 2:T) {
    Y <- sample(1:n,1)
    alpha <- dbinom(Y,n,p)/dbinom(X[t-1],n,p)
    if (runif(1) < alpha){
      X[t] <- Y
    }else{
      X[t] <- X[t - 1]}
  }
  return(X)
}

# calling M-H algorithm and plotting result
test <- metropolisAlgorithm(40,0.5,5000)
par(mfrow=c(2,1))
plot(test, type = "l")
hist(test) breaks = 40)

Typical output (which makes perfect sense):

Problem

My Metropolis-Hastings problem has a stationary binomial distribution, and all proposal distributions q(i,j) are 0.5. With reference to the plot and histogram, should the algorithm be so clearly centered around 0.5, the probability from the binomial distribution? ``` pi <- function(x, n, p){ # returning value from binomial distribution result <- (factorial(n) / (factorial(x) * factorial(n - x))) * p^x * (1 - p)^(n - x) return(result) } metropolisAlgorithm <- function(n, p, T){ # implementation of the algorithm # @n,p binomial parameters # @T number of time steps X <- rep(runif(1),T) for (t in 2:T) { Y <- runif(1) alpha <- pi(X[t - 1], n, p) / pi(Y, n, p) if (runif(1) < alpha) X[t] <- Y else X[t] < X[t - 1] } return(X) } # calling M-H algorithm and plotting result test <- metropolisAlgorithm(40,0.5,5000) par(mfrow=c(2,1)) plot(test, type = "l") hist(test, breaks = 40) ```

Original source