Erik L - 1 year ago 136

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

Recommended for you: Get network issues from **WhatsUp Gold**. **Not end users.**

Answer Source

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

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