Logo Questions Linux Laravel Mysql Ubuntu Git Menu
 

How can I save a JAGS model object in R?

Tags:

r

jags

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.

like image 252
sieste Avatar asked Aug 21 '14 11:08

sieste


1 Answers

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.

like image 150
Woosah Avatar answered Oct 08 '22 09:10

Woosah