r - R : correct results? 中的 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)

典型输出(非常合理):