在 R 中操作 mcmc.list 对象答案

【问题标题】：Manipulating mcmc.list object in R在 R 中操作 mcmc.list 对象
【发布时间】：2016-02-16 17:46:58
【问题描述】：

我使用通过 rjags 调用的 JAGS 来生成 mcmc.list 对象 foldD_samples，其中包含大量随机节点（>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"

干杯，雅各布

【问题讨论】：

rjags 中可能有方法，但我似乎记得你可以做do.call(rbind.data.frame, foldD_samples)。或许胖子使用效率更高data.table::rbindlist
ps 您可以在列表中应用code.samples，而无需强制转换为数据框
感谢 user20650 ！ do.call(rbind.data.frame, foldD_samples) 运作良好。如果这样发布，我很乐意接受这个作为答案。 data.table::rbindlist 不接受 mcmc.list 对象作为输入。还要注意后记中“coda”的假定拼写错误“代码”。
rbindlist(lapply(foldD_samples,as.data.frame)) 可能会工作......如果效率对你很重要。

标签： r data-manipulation mcmc jags

【解决方案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)

【讨论】：

【解决方案2】：

user20650 给出的答案肯定会起作用，但它可能会很慢。另请注意，在撰写本文时，不推荐使用 rbind_list() 来代替 bind_rows()。

我为自己的目的编写的内容将 mcmc.list 转换为“长格式”data.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 个链的 mcmc.list 对象，总共 50,000 行， 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 秒的经过时间。

【讨论】：