在R中使用MCMC Metropolis-Hastings算法对多维后验分布进行采样

Question

我对使用基于Metropolis-Hastings算法的MCMC技术采样后验分布(因此是贝叶斯方法)非常陌生。为此，我使用了R中的mcmc库。我的分布是多维的。为了检查这个metro算法是否适用于多变量分布，我在一个多维的学生t分布(包mvtnorm，函数dmvt)上成功地进行了测试。现在，我想将相同的东西应用于我的多变量分布(2个变量x和y)，但它不起作用；我得到一个错误: error in X，1: Error of dimensions

下面是我的代码：

library(mcmc)
library(mvtnorm)
my.seed <- 123

logprior<-function(X,...)
{
      ifelse( (-50.0 <= X[,1] & X[,1]<=50.0) & (-50.0 <= X[,2] & X[,2]<=50.0), return(0), return(-Inf))
}

logpost<-function(X,...)
{
      log.like <- log( exp(-((X[,1]^2 + X[,2]^2 - 4)/10 )^2) * sin(4*atan(X[,2]/X[,1])) )
      log.prior<-logprior(X)
      log.post<-log.like + log.prior # if flat prior, the posterior distribution is the likelihood one
      return (log.post)
}

x <- seq(-5,5,0.15)
y <- seq(-5,5,0.15)
X<-cbind(x,y)

#out <- metrop(function(X) dmvt(X, df=3, log=TRUE), 0, blen=100, nbatch=100) ; this works
out <- metrop(function(X) logpost(X), c(0,0), blen=100, nbatch=100)
out <- metrop(out)
out$accept 

因此，我尝试尊重与MWE相同的格式，但它仍然不能工作，因为我得到了前面提到的错误。另一件事是，将logpost应用到X可以完美地工作。

提前感谢你的帮助，最好的

Answer 1

metrop函数传递单个样本，因此向logpost传递一个简单的向量，而不是传递一个矩阵(这就是X )。因此，解决方案是将X[,1]和X[,2]分别更改为X[1]和X[2]。

我这样运行它，它导致了其他问题(对于初始化，X[2]/X[1]是NaN )，但这更多地与您特定的可能性模型有关，并且超出了您的问题的范围。