# 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)
## Not run:
# # 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")
# ## End(Not run)
# 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)))
## Not run:
# # 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')
# ## End(Not run)
# Run the same model using 8 chains in parallel:
# distributed computing cluster:
## Not run:
#
# # 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')
# ## End(Not run)
Run the code above in your browser using DataLab