在R中操作mcmc.list对象

Question

我使用了JAGS并通过rjags调用它来生成mcmc.list对象foldD_samples，该对象包含了大量随机节点的trace监视器（>800个节点）。

现在我想使用R来计算这些节点的一个相当复杂的标量值函数，并将输出写入一个mcmc对象，以便我可以使用coda来总结后验分布并运行收敛诊断。

然而，我一直无法弄清楚如何从foldD_samples中获取后验抽样数据到一个数据框中。非常感谢任何帮助。

以下是mcmc.list的结构：

str(foldD_samples)
List of 3
 $ : mcmc [1:5000, 1:821] -0.667 -0.197 -0.302 -0.204 -0.394 ...
  ..- attr(*, "dimnames")=List of 2
  .. ..$ : NULL
  .. ..$ : chr [1:821] "beta0" "beta1" "beta2" "dtau" ...
  ..- attr(*, "mcpar")= num [1:3] 4100 504000 100
 $ : mcmc [1:5000, 1:821] -0.686 -0.385 -0.53 -0.457 -0.519 ...
  ..- attr(*, "dimnames")=List of 2
  .. ..$ : NULL
  .. ..$ : chr [1:821] "beta0" "beta1" "beta2" "dtau" ...
  ..- attr(*, "mcpar")= num [1:3] 4100 504000 100
 $ : mcmc [1:5000, 1:821] -0.492 -0.679 -0.299 -0.429 -0.421 ...
  ..- attr(*, "dimnames")=List of 2
  .. ..$ : NULL
  .. ..$ : chr [1:821] "beta0" "beta1" "beta2" "dtau" ...
  ..- attr(*, "mcpar")= num [1:3] 4100 504000 100
 - attr(*, "class")= chr "mcmc.list"

祝福， Jacob

- Jacob Socolar

---

在rjags中可能有方法，但我记得你可以使用do.call(rbind.data.frame, foldD_samples)。也许更快更高效的方法是使用data.table::rbindlist。 - user20650

请注意，您可以在列表上应用code.samples，而无需将其强制转换为数据框。 - user20650

谢谢user20650！do.call(rbind.data.frame, foldD_samples) 运行得很好。如果以这样的形式发布，我很乐意接受它作为答案。data.table::rbindlist 不接受mcmc.list对象作为输入。还请注意附言中"code"被误写成"coda"的可能打字错误。 - Jacob Socolar

如果效率很重要的话，rbindlist(lapply(foldD_samples,as.data.frame)) 可能会起作用。 - Ben Bolker

Answer 1

由于这是一个 list 结构，您可以使用以下任何一种方法将矩阵绑定在一起。

do.call(rbind.data.frame, foldD_samples)

或者

rbindlist(lapply(foldD_samples, as.data.frame)) # thanks to BenBolker

---

一个示例

library(rjags)
library(coda)
library(data.table)

mod <- textConnection("model {
  A ~ dnorm(0, 1)
  B ~ dnorm(0, 1)
}")

# evaluate
mod <- jags.model(mod, n.chains = 4, n.adapt = 50000) 
pos <- coda.samples(mod,  c("A", "B"),  10000)

out <- do.call(rbind.data.frame, pos)
out2 <- rbindlist(lapply(pos, as.data.frame))
all.equal(out, out2, check.attributes=FALSE)

- user20650

Answer 2

用户20650提供的答案肯定是可行的，但速度可能会很慢。此外，请注意，在撰写本文时，rbind_list（）已被弃用，取而代之的是bind_rows（）。

我为自己的目的编写的某些内容将mcmc.list转换为“长格式”数据.frame。在我的机器上，它比上述方法快4-7倍，并添加了两个额外的列：一个用于链编号，另一个用于步骤编号。

parameter_names <- varnames(mcmc_list)
saved_steps <- as.integer(row.names(mcmc_list[[1]]))
out <- data.frame("chain" = factor(rep(1 : length(mcmc_list), each = length(saved_steps))),
                  "step" = rep(saved_steps, length(mcmc_list)) )
for (param in parameter_names) {
    out[param] <- NA
}
for (a_chain in 1 : length(mcmc_list)) {
    out[out$chain == a_chain, parameter_names ] <- as.data.frame(mcmc_list[[a_chain]])
}
return(out)

使用包含3个链路和50000行的mcmc.list对象， rbind_list/do.call方法：平均0.71秒的经过时间 我的方法：平均0.12秒的经过时间

编辑：对库“coda”中代码的进一步阅读显示，as.matrix()更快。

parameter_names <- varnames(mcmc_list)
saved_steps <- as.integer(row.names(mcmc_list[[1]]))
out <- data.frame("chain" = factor(rep(1 : length(mcmc_list), each = length(saved_steps))),
                  "step" = rep(saved_steps, length(mcmc_list)) )
out <- cbind(out, as.data.frame(as.matrix(chain_samples)))

平均需要0.03秒的时间。

- prince_of_pears