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

Posted 2023-02-19

【中文标题】使用 R 中的 MCMC Metropolis-Hastings 算法对多维后验分布进行采样

【英文标题】：sampling a multimensional posterior distribution using MCMC Metropolis-Hastings algo in R

【发布时间】：2016-10-23 03:52:36

【问题描述】：

我在使用基于 Metropolis-Hastings 算法的 MCMC 技术对后验分布（因此是贝叶斯方法）进行采样方面非常陌生。 为此，我正在使用 R 中的 mcmc 库。我的分布是多维的。为了检查这个 Metro 算法是否适用于多元分布，我在多维学生 T 分布（包 mvtnorm，函数 dmvt）上成功完成了它。 现在我想将相同的东西应用于我的多元分布（2 vars x 和 y），但它不起作用；我收到一个错误：X[, 1] 中的错误：维数不正确

这是我的代码：

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 效果很好。

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

【参考方案1】：

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

我这样运行它，它会导致其他问题（X[2]/X[1] 是初始化的 NaN），但这与您的特定可能性模型有关，超出了您的问题范围。