I am using the package rjags
to do MCMC in R and I would like to save the output of the function jags.model
for later use in another R session.
Here is a simple example for the mean of a Normal distribution:
library(rjags)
N <- 1000
x <- rnorm(N, 0, 5)
model.str <- 'model {for (i in 1:N) {
x[i] ~ dnorm(mu, 5)}
mu ~ dnorm(0, .0001)}'
jags <- jags.model(textConnection(model.str), data = list(x = x, N = N))
update(jags, 1000)
I can generate samples of mu
like this:
coda.samples(model=jags,n.iter=1,variable.names="mu")
# [[1]]
# Markov Chain Monte Carlo (MCMC) output:
# Start = 2001
# End = 2001
# Thinning interval = 1
# mu
# [1,] 0.2312028
#
# attr(,"class")
# [1] "mcmc.list"
Now I would like to save the model object jags
for later use in a new R session, so that I don't have to initialize and burn in the Markov Chain again:
save(file="/tmp/jags.Rdata", list="jags")
quit()
However, after starting a new R session and reloading the model I get an error message that the JAGS model must be recompiled:
load("/tmp/jags.Rdata")
coda.samples(model=jags,n.iter=1,variable.names="mu")
# Error in model$iter() : JAGS model must be recompiled
Why is that? How can I save the object jags
in R for later use?
Note: The question has been asked before, but the OP was not very specific about the problem.
Maybe I am totally off track regarding what you really want to do, but I would set up a jags model like this, using R2jags instead of rjags (just something like a different wrapper):
library(R2jags)
N <- 1000
x <- rnorm(N, 0, 5)
sink("test.txt")
cat("
model{
for (i in 1:N) {
x[i] ~ dnorm(mu, 5)
}
mu ~ dnorm(0, .0001)
}
",fill = TRUE)
sink()
inits <- function() {
list(
mu = dnorm(1, 0, 0.01))
}
params <- c("mu")
chains <- 3
iter <- 1000
jags1 <- jags(model.file = "test.txt", data = list(x = x, N = N),
parameters.to.save = params, inits = inits,
n.chains = chains, n.iter = iter, n.burnin=floor(iter/2),
n.thin = ifelse(floor(iter/100) < 1, 1, floor(iter/100)))
jags2 <- update(jags1, 10000)
jags2
plot(jags2)
traceplot(jags2)
jags2.mcmc <- as.mcmc(jags2)
There is no difference in the results and I like this procedure because it's much more the way I used winbugs, so...
The last line of code converts the jags2-object to an mcmc-list which can be treated by package coda.
Good luck!
P.S. Here's a second answer:
After looking again on your code, the only thing after loading the jags-object that is missing to get the behavior you want, is:
jags$recompile()
coda.samples(model=jags,n.iter=1,variable.names="mu")
But if you really just want to use the already obtained posterior samples or maybe just want to update the chains for more iterations, you can also use the R2jags-procedure.
If you love us? You can donate to us via Paypal or buy me a coffee so we can maintain and grow! Thank you!
Donate Us With