I am using rjags as a sampler. The model has 3 matrices defined. The coda.samples function returns a list of samples. If I take the first sample list the column names look something like this:
> 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]"
Obviously, A, B and C are matrices in my model. I want to reconstruct them based on the mean of these samples. I can easily get the means with colMeans(output[[1]]) but I have no idea how to easily reconstruct the matrices from this vector.
A good way for reconstruction would be the relist() function. So If I had matrices A, B and C in a list L = list(A=A,B=B,C=C) then I could transform this list to a vector with unlist() and convert back with relist(). I am looking for something similar/readymade for mcmc objects, but without avail so far - I can't believe I am the first one needing this. Obviously, relist(colMeans(output[[1]])) doesn't work.
Anyone can help me with reconstructing?
Edit: note also that the relist() function only needs a skeleton, so extracting the skeleton from colnames(output[[1]]) would also do the trick. Or am I complicating?
I don't think relist() will do the trick...
I assume your object output is an object of class mcmc.list, as defined in the R package coda, and output[[1]] is an object of class mcmc representing the sample for the first MCMC chain.
I'm pretty sure coda does not have any understanding that e.g. "A[1,1]" is a JAGS matrix, it just handles it as a variable name. So you'll have to iterate over the relevant variables and impose the structure yourself.
Ideally you'd wrap this in a function like the following:
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]
)
)
)
}
Thus, for example, if your matrix B stored in output[[1]] has 3 rows and 4 columns you'd write:
getMatrix(output[[1]], "B", 3, 4)
to get the means as a matrix object in R.