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

Question

我很新于使用基于Metropolis-Hastings算法的MCMC技术对后验分布进行抽样（因此是贝叶斯方法）。我在R中使用mcmc库来实现。我的分布是多维度的。为了检查这个metro算法是否适用于多元分布，我成功地在一个多维学生t分布上进行了尝试（使用mvtnorm包中的dmvt函数）。

现在我想将同样的方法应用到我的多变量分布（2个变量x和y），但它不起作用。我得到了一个错误：Error in X[,1]: incorrect number 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完美地运行。

感谢您的帮助，最好的祝福。

- Jean-Philippe Fontaine

Answer 1

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

我按照这种方式运行它，但会导致其他问题（初始化时X [2] / X [1]为NaN），但这与您特定的似然模型有关，超出了您的问题范围。

- merv