r - 需要有关 R 中 MCMC 估计的建议

Question

假设随机变量 X ∼ N(μ,σ2) 分布
我得到了观察结果 x =(0.7669, 1.6709, 0.8944, 1.0321, 0.0793, 0.1033, 1.2709, 0.7798, 0.6483, 0.3256)
给定先验分布 μ ∼ N(0, (100)^2) 和 σ2 ∼ Inverse − Gamma(5/2, 10/2)。
我正在尝试对参数进行 MCMC 估计。
这是我的示例代码，我尝试使用 Random Walk MCMC 进行估计，但似乎没有收敛:

x = c(0.7669,1.6709,0.8944, 1.0321, 0.0793, 0.1033, 1.2709, 0.7798, 0.6483, 0.3256)
mu.old = rnorm(1,0,100); #initial guess for mu
alpha = 5/2;beta = 10/2
sigma.old = sqrt(rinvgamma(1, shape = alpha, scale = 1/beta)) #initial guess for sigma
burn.in = 50000
rep = 100000
acc = 0

lik.old = sum(dnorm(x,mu.old,sd = sigma.old),log = T)
res = c(mu.old,sigma.old,lik.old)
for (i in (1:rep)){
  mu.new = mu.old + rnorm(1,0,10)
  sigma.new = sigma.old + rnorm(1,0,10)
  if ((mu.new <=0)||sigma.new <= 0){
    res = rbind(res,c(mu.old,sigma.old,lik.old))
  }
  else{
    lik.new = sum(dnorm(x,mu.new,sigma.new,log = T))
    r = exp(lik.new - lik.old)
    u = runif(1,0,1)
    if (u<r){
      res = rbind(res,c(mu.new,sigma.new,lik.new))
      mu.old = mu.new
      sigma.old = sigma.new
      lik.old = lik.new
      acc = acc + 1
    }
    else{
      res = rbind(res,c(mu.old,sigma.old,lik.old))
    }
  }
  if (i %% 10000 == 0){
    plot(res[,1],type = 'l')
  }
}

有人可以就如何找到更好的 MCMC 方法提出建议吗？

Answer 1

代码有几个问题。

您计算对数似然，但忽略对数先验。如此有效地您正在使用统一的先验，而不是您指定的那个。

对数似然的初始计算在错误的地方有括号:

sum(dnorm(x,mu.old,sd = sigma.old),log = T)

最后一个参数应该是 log = TRUE ，内dnorm()打电话，不在外面。 (避免将 T 用于 TRUE 。键入时节省 3 个字符是不值得的。)

使用 rbind()发展您的 res阵列很慢。最好提前分配整个事情，并将结果分配给特定的行，例如

# initialize
res <- matrix(NA, nrow = reps, ncol = 3)

# in the loop
res[i,] <- c(mu.old,sigma.old,lik.old)