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)

典型输出（非常合理）：

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

Metropolis-Hastings algorithm in R: correct results?

r

markov-chains