snowで並列化してJAGSを実行する関数 [統計]
自作のJAGS実行R関数をsnow対応にしてみた。例によってエラーチェックとかはいいかげん。あと、JAGS 2対応は不十分。
使用例
library(coda) library(snow) 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 & is.vector(dat)) { 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 if (is.matrix(dat) || is.array(dat)) { dat.l <- paste(dat.l, "structure(c(") v <- as.vector(dat) for (j in 1:(length(v) - 1)) { dat.l <- paste(dat.l, v[j], ", ", sep = "") } dat.l <- paste(dat.l, dat[length(v)], "),", sep = "") dat.l <- paste(dat.l, ".Dim=c(", paste(dim(dat), collapse = ","), "))") } else { if (is.character(dat)) { dat <- paste("\"", dat, "\"", sep = "") } 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, (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) if (n.chains == 1) { write(paste("parameters in \"", init.file, "\", chain(", 1, ")", sep=""), cmd.file, append = TRUE) } else { 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", format(n.burnin, scientfic = FALSE)), cmd.file, append = TRUE) for (i in 1:length(parameters.to.save)) { write(paste("monitor set ", parameters.to.save[i], ",thin(", format(n.thin, scientific = FALSE), ")", sep=""), cmd.file, append = TRUE) } write(paste("update", format(n.iter - n.burnin, scientific = FALSE)), 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, (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 larger than 0.") } else if (n.chains >= 2 & n.chains != length(inits)) { stop("n.chains != length(inits).") } else { result <- vector("list", n.chains) jags.makedata(data = data, file = data.file) if (n.chains == 1) { jags.makeinit(inits = inits, file = init.file) } else { 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) if (system(paste(jags, cmd.file, sep = " "), wait = TRUE) == 0) { for (i in 1:n.chains) { output.file <- paste(out, "chain", i, ".txt", sep = "") index.file <- paste(out, "index.txt", sep = "") result[[i]] <- read.coda(output.file, index.file) } rslt <- mcmc.list(result) } else { stop("calculation failed!") } } return(rslt) } jags.parrun <- function(data, inits, parameters.to.save, model.file = "jags.bug", data.prefix = "jags-data", init.prefix = "jags-init", cmd.prefix = "jags", out = "CODA", jags = "/usr/local/bin/jags", n.chains = 3, n.iter = 2000, n.burnin = floor(n.iter/2), n.thin = max(1, (n.iter - n.burnin) %/% 1000), spec = rep("localhost", n.chains)) { if (!file.exists(model.file)) { stop(paste(model.file, "does not exist.")) } if (n.chains < 1) { stop("n.chains must be larger than 0.") } else if (n.chains >= 2 & n.chains != length(inits)) { stop("n.chains != length(inits).") } else if (!is.list(data)) { stop("data must be a list.") } else { par <- parameters.to.save iter <- n.iter burnin <- n.burnin thin <- n.thin params <- list("jags.run", "mcmc", "mcmc.list", "read.coda", "jags.makedata", "jags.makeinit", "jags.makecmd") cl <- makeCluster(spec, type = "SOCK") clusterExport(cl, params) r <- parLapply(cl, 1:n.chains, function(i) { jags.run(data, inits[[i]], par, model.file = model.file, data.file = paste(data.prefix, i, ".R", sep = ""), init.file = paste(init.prefix, i, ".R", sep = ""), cmd.file = paste(cmd.prefix, i, ".R", sep = ""), out = paste(out, i, sep = ""), jags = jags, n.chains = 1, n.iter = iter, n.burnin = burnin, n.thin = thin) }) stopCluster(cl) rslt <- mcmc.list(lapply(1:n.chains, function(x) r[[x]][[1]])) } return(rslt) }
使用例
read.jagsdata <- function(file) { e <- new.env() eval(parse(file), e) return(as.list(e)) } ## set data file <- "blocker.bug" dt <- read.jagsdata("blocker-data.R") ## set inits ini <- list() ini[[1]] <- ini[[2]] <- ini[[3]] <- read.jagsdata("blocker-init.R") ## set RNGs ini[[1]]$.RNG.name <- "base::Mersenne-Twister" ini[[1]]$.RNG.seed <- 123 ini[[2]]$.RNG.name <- "base::Mersenne-Twister" ini[[2]]$.RNG.seed <- 1234 ini[[3]]$.RNG.name <- "base::Mersenne-Twister" ini[[3]]$.RNG.seed <- 12345 ## set parameters vars <- c("d", "delta.new", "sigma") iter <- 30000 adapt <- 3000 thin <- 10 chains <- 3 post <- jags.parrun(data = dt, inits = ini, parameters.to.save = vars, model.file = file, n.chains = chains, n.iter = iter + adapt, n.burnin = adapt, n.thin = thin)
コメント 0