我通常从具有多个链的 rjags 中调用 JAGS 用于诊断目的(例如,4 个链)。之后我经常想对后验参数估计做一些后处理(例如,使用预测值,计算额外的统计数据)。但是,此时将链存储在列表中是一件麻烦事。
将链组合成单个参数列表的好方法是什么?
包runjags
有功能combine.mcmc
。它的默认设置是组合一个或多个链并返回单个链。例如,
library(runjags)
fit <- combine.mcmc(multichainfit)
它还有其他组合链的选项。
这是rjags
和coda
解决方案。假设您生成两个这样的链:
library(rjags)
x <- rnorm(100)
model.str <- 'model {for (i in 1:100) {
x[i] ~ dnorm(mu, sd)}
mu ~ dnorm(0, .0001)
sd ~ dgamma(0.1, 0.1)}'
jags <- jags.model(textConnection(model.str), data = list(x = x),n.chains=2)
smpls <- coda.samples(model=jags,n.iter=2,variable.names=c("mu","sd"))
smpls
# [[1]]
# Markov Chain Monte Carlo (MCMC) output:
# Start = 1
# End = 2
# Thinning interval = 1
# mu sd
# [1,] -0.09152588 0.9009973
# [2,] 0.05586651 1.0482184
#
# [[2]]
# Markov Chain Monte Carlo (MCMC) output:
# Start = 1
# End = 2
# Thinning interval = 1
# mu sd
# [1,] 0.06893182 0.7317955
# [2,] 0.13599206 0.8517304
#
# attr(,"class")
# [1] "mcmc.list"
您可以通过以下方式将两个(或任意多个)链组合成一个矩阵
do.call(rbind, smpls)
# mu sd
# [1,] -0.09152588 0.9009973
# [2,] 0.05586651 1.0482184
# [3,] 0.06893182 0.7317955
# [4,] 0.13599206 0.8517304
如果你想要一个类的对象mcmc
,只需使用函数mcmc
:
mcmc(do.call(rbind, smpls))
# [[1]]
# Markov Chain Monte Carlo (MCMC) output:
# Start = 1
# End = 4
# Thinning interval = 1
# mu sd
# [1,] -0.09152588 0.9009973
# [2,] 0.05586651 1.0482184
# [3,] 0.06893182 0.7317955
# [4,] 0.13599206 0.8517304
如果需要,您可以保留原始链的start
、end
或thin
属性。例如保存thin
:
mcmc(do.call(rbind, smpls), thin=thin(smpls))
(当然你不能保留所有三个属性)