SSブログ

JAGS 0.99.0用Rコード [統計]

JAGS 0.99.0をRから使う関数。例によってエラー処理などは手を抜きまくりである。本家本元のrjagsがJAGS 1.0.0リリース後に公開される予定とのことなので、その後はそちらを利用することになるだろう。

library(coda)

jags.makedata <- function(data, file = "jags-data.R") {
  if (is.list(data)) {
    for (i in 1:length(data)) {
      lab <- labels(data)[i]
      dat.l <- paste("\"", lab, "\" <- ", sep="")
      dat <- data[[i]]
      if (length(dat) > 1) {
        dat.l <- paste(dat.l, "c(")
        for (j in 1:(length(dat) - 1)) {
          dat.l <- paste(dat.l, dat[j], ", ", sep="")
        }
        dat.l <- paste(dat.l, dat[length(dat)], ")", sep="")
      } else {
        dat.l <- paste(dat.l, as.character(dat))
      }
      write(dat.l, file, append=(i > 1))
    }
  } else {
    stop("jags.makedata: first parameter must be a list.")
  }
}

jags.makeinit <- function(inits, file = "jags-init.R") {
  jags.makedata(inits, file)
}

jags.makecmd <- function(parameters.to.save,
                         cmd.file = "jags.cmd",
                         model.file = "jags.bug",
                         data.file = "jags-data.R",
                         init.file = "jags-init.R",
                         out = "CODA",
                         n.chains = 3,
                         n.iter = 2000,
                         n.burnin = floor(n.iter/2),
                         n.thin = max(1, floor((n.iter - n.burnin)/1000))) {
  write(paste("model in \"", model.file, "\"", sep=""), cmd.file)
  write(paste("data in \"", data.file, "\"", sep=""), cmd.file,
        append = TRUE)
  write(paste("compile, nchains(", n.chains, ")", sep=""), cmd.file,
        append = TRUE)
  for (i in 1:n.chains) {
    init.i <- gsub("^([^\\.]*)\\.([^\\.]*$)",
                   paste("\\1-", i, ".\\2", sep = ""),
                   init.file)
    write(paste("parameters in \"", init.i, "\", chain(", i, ")",
                sep=""), cmd.file,
          append = TRUE)
  }
  write("initialize", cmd.file, append = TRUE)
  write(paste("update", as.character(n.burnin)), cmd.file,
        append = TRUE)
  for (i in 1:length(parameters.to.save)) {
    write(paste("monitor set ", parameters.to.save[i],
                ",thin(", as.character(n.thin), ")", sep=""), cmd.file,
          append = TRUE)	
  }
  write(paste("update", as.character(n.iter - n.burnin)), cmd.file,
        append = TRUE)
  write(paste("coda *,stem(", out, ")", sep=""), cmd.file,
        append = TRUE)
  write("exit", cmd.file, append = TRUE)
}

jags.run <- function(data, inits, parameters.to.save,
                     model.file = "jags.bug",
                     data.file = "jags-data.R",
                     init.file = "jags-init.R",
                     cmd.file = "jags.cmd",
                     out = "CODA",
                     jags = "/usr/local/bin/jags",
                     n.chains = 3,
                     n.iter = 2000,
                     n.burnin = floor(n.iter/2),
                     n.thin = max(1, floor((n.iter - n.burnin)/1000))) {
  if (!file.exists(model.file)) {
    stop(paste(model.file, "does not exist."))
  }
  if (n.chains < 1) {
  	stop("n.chains must be equal to or larger than 0.")
  } else  {
    result <- vector("list", n.chains)
    jags.makedata(data = data, file = data.file)
    for (i in 1:n.chains) {
      init.i <- gsub("^([^\\.]*)\\.([^\\.]*$)",
                     paste("\\1-", i, ".\\2", sep = ""),
                     init.file)
      jags.makeinit(inits = inits[[i]], file = init.i)
    }
    jags.makecmd(parameters.to.save = parameters.to.save,
                 cmd.file = cmd.file,
                 model.file = model.file,
                 data.file = data.file,
                 init.file = init.file,
                 out = out,
                 n.chains = n.chains,
                 n.iter = n.iter,
                 n.burnin = n.burnin,
                 n.thin = n.thin)
    system(paste(jags, cmd.file, sep = " "), wait = TRUE)
    for (i in 1:n.chains) {
      result[[i]] <- read.coda(output.file = paste(out, "chain", i, ".txt", sep=""),
                               index.file = paste(out, "index.txt", sep = ""))
    }
    rslt <- mcmc.list(result)
  }
  return(rslt)
}
実行例
model.file <- system.file(package = "R2WinBUGS", "model", "schools.txt")
data(schools)
J <- nrow(schools)
y <- schools$estimate
sigma.y <- schools$sd
data <- list (J = J, y = y, sigma.y = sigma.y)
init1 <- list(.RNG.name = "\"base::Mersenne-Twister\"",
    .RNG.seed = 11,
    theta = rnorm(J, 0, 100),
    mu.theta = rnorm(1, 0, 100),
    sigma.theta = runif(1, 0, 100))
init2 <- list(.RNG.name = "\"base::Mersenne-Twister\"",
    .RNG.seed = 113,
    theta = rnorm(J, 0, 100),
    mu.theta = rnorm(1, 0, 100),
    sigma.theta = runif(1, 0, 100))
init3 <- list(.RNG.name = "\"base::Mersenne-Twister\"",
    .RNG.seed = 1131,
    theta = rnorm(J, 0, 100),
    mu.theta = rnorm(1, 0, 100),
    sigma.theta = runif(1, 0, 100))
parameters <- c("theta", "mu.theta", "sigma.theta")

schools.sim <- jags.run(data,
    inits = list(init1, init2, init3),
 	parameters.to.save = parameters,
    model = model.file,
    n.chains = 3, n.iter = 5000)

nice!(0)  コメント(0)  トラックバック(1) 
共通テーマ:パソコン・インターネット

nice! 0

コメント 0

コメントを書く

お名前:
URL:
コメント:
画像認証:
下の画像に表示されている文字を入力してください。

Facebook コメント

トラックバック 1