runjags.options(new.windows=FALSE)
# run a model to calculate the intercept and slope of the expression
# y = m x + c, assuming normal observation errors for y:
# Simulate the data
X <- 1:100
Y <- rnorm(length(X), 2*X + 10, 1)
# Model in the JAGS format
model <- "model {
for(i in 1 : N){
Y[i] ~ dnorm(true.y[i], precision);
true.y[i] <- (m * X[i]) + c
}
m ~ dunif(-1000,1000)
c ~ dunif(-1000,1000)
precision ~ dexp(1)
}"
# Data and initial values in a named list format,
# with explicit control over the random number
# generator used for each chain (optional):
data <- list(X=X, Y=Y, N=length(X))
inits1 <- list(m=1, c=1, precision=1,
.RNG.name="base::Super-Duper", .RNG.seed=1)
inits2 <- list(m=0.1, c=10, precision=1,
.RNG.name="base::Wichmann-Hill", .RNG.seed=2)
# Run the model and produce plots
results <- run.jags(model=model, monitor=c("m", "c", "precision"),
data=data, n.chains=2, method="rjags", inits=list(inits1,inits2))
# Standard plots of the monitored variables:
plot(results)
# Look at the summary statistics:
print(results)
# Extract only the coefficient as an mcmc.list object:
library('coda')
coeff <- as.mcmc.list(results, vars="m")
# The same model but using embedded shortcuts to specify data, inits and monitors,
# and using parallel chains:
# Model in the JAGS format
model <- "model {
for(i in 1 : N){ #data# N
Y[i] ~ dnorm(true.y[i], precision) #data# Y
true.y[i] <- (m * X[i]) + c #data# X
}
m ~ dunif(-1000,1000) #inits# m
c ~ dunif(-1000,1000)
precision ~ dexp(1)
#monitor# m, c, precision
}"
# Simulate the data
X <- 1:100
Y <- rnorm(length(X), 2*X + 10, 1)
N <- length(X)
initfunction <- function(chain) return(switch(chain,
"1"=list(m=-10), "2"=list(m=10)))
# Run the 2 chains in parallel (allowing the run.jags function
# to control the number of parallel chains). We also use a
# mutate function to convert the precision to standard deviation:
results <- run.jags(model, n.chains=2, inits=initfunction,
method="parallel", mutate=list("prec2sd", vars="precision"))
# View the results using the standard print method:
results
# Look at some plots of the intercept and slope on a 3x3 grid:
plot(results, c('trace','histogram','ecdf','crosscorr','key'),
vars=c("m","^c"),layout=c(3,3))
# Write the current model representation to file:
write.jagsfile(results, file='mymod.txt')
# And re-run the simulation from this point:
newresults <- run.jags('mymod.txt')
# Run the same model using 8 chains in parallel:
# distributed computing cluster:
# A list of 8 randomly generated starting values for m:
initlist <- replicate(8,list(m=runif(1,-20,20)),simplify=FALSE)
# Run the chains in parallel using JAGS (2 models
# with 4 chains each):
results <- run.jags(model, n.chains=8, inits=initlist,
method="parallel", n.sims=2)
# Set up a distributed computing cluster with 2 nodes:
library(parallel)
cl <- makeCluster(4)
# Run the chains in parallel rjags models (4 models
# with 2 chains each) on this cluster:
results <- run.jags(model, n.chains=8, inits=initlist,
method="rjparallel", cl=cl)
stopCluster(cl)
# For more examples see the quick-start vignette:
vignette('quickjags', package='runjags')
# And for more details about possible methods see:
vignette('userguide', package='runjags')Run the code above in your browser using DataLab