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)

that's a good question, I did not even consider to look for the native function :/ Thanks for the tip! — G. Debailly, Oct 30 '16 at 13:00
When I repeatedly run your code I typically don't see what you claim. I see something that seems to be concentrated in the vicinity of the first random number generated, which is often quite far from 0.5. Thus, there is a problem with your implementation. Your proposal distribution seems odd. — John Coleman, Oct 30 '16 at 13:17
Shouldn't `X` take on integer values rather than values in `[0,1]`? — John Coleman, Oct 30 '16 at 13:38
`X[t] < X[t - 1]` should be `X[t] <- X[t - 1]`, though there is the semantic issue that `X` seems to be the number of successes (discrete in 1:40) but you are trying to simulate it via a random walk in `[0,1]` — John Coleman, Oct 30 '16 at 14:14
dear all, thanks so much for attempting to help. It seems my problem is really a lack of understanding the problem. And @JohnColeman, thanks for both noticing the typo and addressing what seem to be the fundamental problem. — G. Debailly, Oct 30 '16 at 14:23
meanwhile, I would like to up vote all your comments. I think I'm two rep points short of doing that :) — G. Debailly, Oct 30 '16 at 14:26
guys, I'm sorry to admit that I properly mistook my proposal distribution. That distribution actually needed its own function, returning Y based on three different states of X(t-1). The solution was sitting right in front of my nose, I just did not understand it. Thanks again! — G. Debailly, Oct 30 '16 at 15:56

score 4 · Accepted Answer · answered Oct 30 '16 at 14:28

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

you actually understood the problem even though I did not present it correctly/ completely. What a hero, thank you so much! — G. Debailly, Oct 30 '16 at 14:37
@ErikL By a sheer coincidence, I was sitting down at the computer to work on an implementation of the Metropolis algorithm (in Excel VBA, though I have previously done it in R as well) when I saw this post. — John Coleman, Oct 30 '16 at 14:40
that coincidence made good. I just love this place in all its tail-saving glory :) — G. Debailly, Oct 30 '16 at 15:53

Metropolis-Hastings algorithm in R: correct results?

1 Answers1