从mcmc对象重构变量

Question

我正在使用rjags作为取样器。该模型定义了3个矩阵。coda.samples函数返回示例列表。如果我使用第一个示例列表，列名如下所示：

> colnames(output[[1]])
"A[1,1]"  "A[2,1]"  "A[1,2]"  "A[2,2]" ... 
"B[1,1]"  "B[2,1]"  "B[3,1]"  "B[4,1]" ... 
"C[1,1]"  "C[2,1]"

显然，在我的模型中，A，B和C是矩阵。我想根据这些样本的平均值重建它们。我可以很容易地得到colMeans(output[[1]])的方法，但是我不知道如何从这个向量中很容易地重建矩阵。

relist()函数是一种很好的重建方法。因此，如果在一个列表L = list(A=A,B=B,C=C)中有矩阵A、B和C，那么我可以用unlist()将这个列表转换为向量，然后用relist()进行转换。我正在寻找类似的东西/现成的mcmc对象，但至今没有效果-我不敢相信我是第一个需要这个。显然，relist(colMeans(output[[1]]))不起作用。

有人能帮我重建吗？

编辑:还请注意，relist()函数只需要一个骨架，所以从colnames(output[[1]])中提取骨架也是很有用的。还是我变得复杂了？

Answer 1

我不认为relist()会做这个.

我假设您的对象output是mcmc.list类的对象，如R包coda中定义的那样，而output[[1]]是类mcmc的对象，表示第一个MCMC链的示例。

我确信coda没有任何理解，例如，"A[1,1]"是一个JAGS矩阵，它只是将它作为一个变量来处理。因此，您必须迭代相关变量，并自己强制执行结构。

理想情况下，您应该将其包装在如下函数中：

getMatrix <- function(output, varname, rows, cols) {
  unname(
    sapply(1:cols, function(j)
      sapply(1:rows, function(i)
        summary(output[,sprintf("%s[%s,%s]", varname, i, j)])[[1]][1]
      )
    )
  )
}

因此，例如，如果存储在B中的矩阵output[[1]]有3行4列，那么您应该编写：

getMatrix(output[[1]], "B", 3, 4)

在R中得到作为矩阵对象的方法。