2017-04-11 82 views
0

我使用rjags作爲採樣器。該模型定義了3個矩陣。函數coda.samples返回樣本列表。如果我拿第一個樣本列表,列名看起來像這樣:重建mcmc對象的變量

> 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]])提取骨架也可以做到這一點。還是我複雜?

回答

0

我不認爲relist()會做的伎倆......

我假設你的對象outputmcmc.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] 
    ) 
    ) 
) 
} 

因此,舉例來說,如果你存儲在output[[1]]矩陣B有3行4列,你會寫:

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

獲得作爲R中的矩陣對象的手段。