Erik L Erik L - 1 month ago 11
R Question

Metropolis-Hastings algorithm in R: correct results?

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)


enter image description here

Answer

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):

enter image description here