StatModeling Memorandum

StatModeling Memorandum

StanとRとPythonでベイズ統計モデリングします. たまに書評.

JAGSの結果からbugsオブジェクトを作成する

JAGSを使うとcoda.samples関数で最終的にmcmc.listが得られます。ここからbugsオブジェクトに変換するスクリプトを作りました。

{rjags}, {R2WinBUGS}の2つのライブラリを読んでおく必要があります。

library(rjags)
library(R2WinBUGS)

mcmc.list2bugs <- function(mcmc.list) {
   b1 <- mcmc.list[[1]]
   m1 <- as.matrix(b1)
   mall <- matrix(numeric(0), 0, ncol(m1))
   n.chains <- length(mcmc.list)
   for (i in 1:n.chains) {
      mall <- rbind(mall, as.matrix(mcmc.list[[i]]))
   }
   sims.array <- array(mall, dim = c(nrow(m1), n.chains, ncol(m1)))
   dimnames(sims.array) <- list(NULL, NULL, colnames(m1))
   mcpar <- attr(b1, "mcpar")
   as.bugs.array(
      sims.array = sims.array,
      model.file = NULL,
      program = NULL,
      DIC = FALSE,
      DICOutput = NULL,
      n.iter = mcpar[2],
      n.burnin = mcpar[1] - mcpar[3],
      n.thin = mcpar[3]
   )
}

これでおもむろに変換したオブジェクトをprintすれば、WinBUGSを使ったときのようにパラメータの分布の要約やGelmanとRubinのR hatなどが表示されます。