Erik L - 1 year ago 136
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)
``````

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

Recommended from our users: Dynamic Network Monitoring from WhatsUp Gold from IPSwitch. Free Download