需要有关 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)
```

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

Needing advice on MCMC estimation in R

estimation

r

mcmc