Nothing
### Initial library
library(pbdMPI, quietly = TRUE)
COMM.SIZE <- comm.size()
COMM.RANK <- comm.rank()
### Set configurations.
I.b <- 1000 # burn-in iteration
I.t <- 10 # thinning iteration
I.n <- 1000 # total samples
### Set data and parameters.
x <- MASS::galaxies # obtain galaxies
n <- length(x)
mu.0 <- mean(x) # prior mean of mu.x
sigma.0 <- sqrt(var(x) / n) # prior std of sigma.x
sigma.x <- sqrt(var(x)) # data std
# Load JAGS
options(warn = -1)
suppressMessages(library(rjags, quietly = TRUE))
options(warn = 0)
load.module("lecuyer", quiet = TRUE)
### JAGS model.
model <- "
model{
### likelihood
for(i in 1:n){
x[i] ~ dnorm(mu, tau.x)
}
tau.x <- 1 / sigma.x
### prior
mu ~ dnorm(mu.0, tau.0)
tau.0 <- 1 / sigma.0
}
"
fn.model <- tempfile()
cat(model, file = fn.model)
### JAGS's monitoring variables.
monitor.vars <- "mu"
### JAGS's data.
data <- list(x = x, n = n, mu.0 = mu.0, sigma.0 = sigma.0, sigma.x = sigma.x)
### Use task.pull.
FUN <- function(jid, n.job){
### JAGS's initial values.
set.seed(123 + jid) # seed for initial rnorm only.
inits <- list(list(mu = rnorm(1, mean = mu.0),
.RNG.name = "lecuyer::RngStream",
.RNG.seed = 1234 + jid)) # seed for JAGS.
### Run rjags.
m <- jags.model(fn.model, data, inits, n.chains = length(inits),
quiet = TRUE)
update(m, n.iter = I.b, progress.bar = "none")
I.n.sub <- ceiling(I.n / n.job) * I.t
x.post <- coda.samples(m, monitor.vars, n.iter = I.n.sub, thin = I.t,
progress.bar = "none")
x.post[[1]]
}
### Run tasks.
n.job <- COMM.SIZE
# ret.all <- pbdLapply(1:n.job, FUN, n.job)
ret.all <- task.pull(1:n.job, FUN, n.job)
### Summary.
if(COMM.RANK == 0){
save(ret.all, file = "tmp.rda")
class(ret.all) <- "mcmc.list"
pdf("galaxy_rjags.pdf", width = 8, height = 4)
plot(ret.all)
dev.off()
}
### Finalize jobs.
finalize()
Any scripts or data that you put into this service are public.
Add the following code to your website.
For more information on customizing the embed code, read Embedding Snippets.