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)
コメント 0