用Metropolis-Hastings算法(MCMC)模拟概率模型

Question

我在R中工作，下面是我的代码：

n <- 1000
trueB0 <- 0
trueB1 <- 0.1
prior.mean0 <- 0
prior.sd0 <- 10
prior.mean1 <- 0
prior.sd1 <- 10

x <- rnorm( n, 0, 1 )
eta <- trueB0 + trueB1 * x
P <- pnorm(eta)
Y <- rbinom(n, 1, P)

prior <- function(param){
  b0=param[1]
  b1=param[2]
  B0prior <- dnorm(b0, mean = prior.mean0, sd = prior.sd0, log = T)
  B1prior <- dnorm(b1, mean = prior.mean1, sd = prior.sd1, log = T)
  return(B0prior + B1prior)
}

likelihood <- function(param){
  b0=param[1]
  b1=param[2]
  eta.l <- b0 + b1 * x
  P.l <- pnorm(eta.l)
  singlelikelihoods = dbinom(Y, 1, P.l, log = T)
  sumall=sum(singlelikelihoods)
  return(sumall)
}

posterior <- function(param){
  return(prior(param)+likelihood(param))
}

rho <- 0
sd0 <- 1
sd1 <- 1
sigma <- matrix(c(sd0^2, sd0*sd1*rho, sd0*sd1*rho, sd1^2), ncol = 2)

proposalfunction <- function(param){
  return(rmvnorm(1, mean = param, sigma))
}

q <- function(m1, m2){
  return(dmvnorm(m1, mean = m2, sigma, log = T))
}

run_metropolis_MCMC <- function(initialvalue, iterations){
  chain = array(dim = c(iterations+1,2))
  chain[1,] = initialvalue
  for (i in 1:iterations){
    proposal = proposalfunction(chain[i, ])
    #bug here in the q evaluatoin part
    alpha = exp(posterior(proposal) - posterior(chain[i, ]) + q(chain[i, ], proposal) - q(proposal, chain[i, ]))
    if (runif(1) < min(1,alpha)) 
    {
      chain[i+1,] = proposal
    }
    else
    {
      chain[i+1,] = chain[i,]
    }
  }  
  return(chain)
}

initialvalue <- c(0 , 0.1)
iterations <- 10000
MH = run_metropolis_MCMC(initialvalue, iterations)

我猜问题是，当计算MH部分的α时，我应该得到q( beta | beta* )/q(beta*|beta) (q是建议函数，beta是当前链值，beta*是Q使用beta生成的值)。然而，我真的不知道如何才能在R代码中得到它。(我知道Q应该是一个二元正态分布密度，但我不知道如何在函数中实现beta和beta*并得到Q的值)

Answer 1

我运行你的代码，如果我添加库(Mvtnorm)，我在Q求值中没有得到任何bug。我实际上搞不懂你的问题是什么，因为你的Q函数给出了确切的结果: q(chaini，，proposal) = q(beta|beta*)，对于二元正态分布(它给出了它的对数密度)。

以下是代码的几个一般注释：

• runif(1) < min(1,alpha) 与 runif(1) < alpha 相同
• 如果您的代码计算已接受提案的数量，那就太好了，这样您就知道其中的百分比。 如果该百分比很小，那么您的提案 sd 就很大，反之亦然。 在 for 循环之前添加 n.accept <-0 并在第一个 if 语句中添加 n.accept <- n.accept + 1 。 然后在返回之前添加

print(c("beta 接受抽签的百分比：", 100*n.accept/iterations))

如果此答案对您有帮助，请将其标记为已接受。