R中的Metropolis-Hastings算法：正确的结果？

Question

我的 Metropolis-Hastings 问题有一个平稳的二项分布，并且所有提案分布 q(i,j) 都是 0.5。参考图和直方图，算法是否应该如此清晰地以二项分布的概率0.5为中心？

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

您有 3 个问题：

1) 您似乎想要模拟二项式分布，因此您的随机游走应该在1:n 范围内的整数上，而不是在[0,1] 范围内的实数上。

2) 你在alpha的计算中切换了分子和分母

3) 你在X[t] < X[t - 1] 中有错字。

修复这些问题并稍微清理一下您的代码（包括按照@ZheyuanLi 的建议使用dbinom 函数）产生：

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)

典型输出（非常有意义）：