Metropolis-Hastings algorithm in R: correct results?

Question

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)

Answer 1

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