Run Parallel MCMC Chains with JAGS/rjags on multiple cores?

[Katherine Lockhart]

Hi Everyone,

This is my first time requesting help from the group, apologies if I do not follow standard procedures. 

I am running some Bayesian models with JAGS/rjags. I have a lot of data in my model (2200 rows) and the chains are relatively long. I want to run each chain on a separate core in order to speed up the model runs. Currently each run is taking about 4 hours. I've searched a bit online, but wound up quite confused after reading about setting up RNGs for each chain, etc. Now, I thought I heard somewhere that rjags does this automatically with multiple chains, but when I look at my CPU monitor it appears that only one core is running. 

Here is my R code for the JAGS portion. Please let me know if you have any advice.

data <- list("y" = nit, "X"=X.mod, "n" = n, "l" = l) 

# Exponential model with gamma priors 

jags.exponential.gamma <-(ExponentialGammaPriors.JAG, data=data, n.chains=2)

update(jags.exponential.gamma, n.iter=100000)

samps <- coda.samples(jags.exponential.gamma, 

                      variable.names=c("beta", "mu"), 

                      n.iter=200000, 

                      thin=50)

Thank you!

Katie 

--

--

Katherine Lockhart, MS

PhD Candidate

Hydrologic Sciences Graduate Group

UC Davis

[Yee, Julie]

Hi, Katie,

I've found this example, courtesy of Tim Handley, is a useful step-by-step template for multicore processing of Bayesian models using rjags:  http://sourceforge.net/p/mcmc-jags/discussion/610036/thread/585b0e4c/

Julie

--
Check out our R resources at http://www.noamross.net/davis-r-users-group.html
---
You received this message because you are subscribed to the Google Groups "Davis R Users' Group" group.
To unsubscribe from this group and stop receiving emails from it, send an email to davis-rug+...@googlegroups.com.
Visit this group at http://groups.google.com/group/davis-rug.
For more options, visit https://groups.google.com/d/optout.

--

____________________________________________________

Julie L. Yee, Statistician

U.S. Geological Survey

Western Ecological Research Center

www.werc.usgs.gov

===========================

*** NEW PHONE and FAX ***

Voice  530 - 669 - 5097

Fax     707 - 678 - 5039

===========================

[Carl Boettiger]

Hi Katie, Julie,

I think the example using the snow package looks a bit overkill if you just want to use multiple cores (but handy for clusters of multiple machines using MPI mode instead of SOCK mode shown in that example.)

Perhaps you were thinking of the `jags.parallel` function from the R2jags package? ?jags.parallel has some examples.  

Cheers,

Carl

Carl Boettiger
UC Santa Cruz
http://carlboettiger.info/

[Katherine Lockhart]

Hi Carl and Julie,

Thank you for your feedback! jags.parallel in R2jags seemed pretty straightforward so I went ahead and tried it. I can get the code to run on one core, but no two. Here is my code and error:

# this runs on 1 core

samps <- jags(data=data,parameters.to.save=c("beta", "mu"),n.iter=3000,

                       model.file="D:/Dropbox/PhD/Loading/Bayesian_Analysis/ExponentialGammaPriors.JAG",

                       n.chains=2, n.burnin=1000,n.thin=5)

# I get an error when I run this one in parallel

samps <- jags.parallel(data=data,parameters.to.save=c("beta", "mu"),n.iter=3000,

              model.file="D:/Dropbox/PhD/Loading/Bayesian_Analysis/ExponentialGammaPriors.JAG",

              n.chains=2, n.burnin=1000,n.thin=5)

Error in get(name, envir = envir) : invalid first argument

#I get the same error if I write my model as a function and then use model.file=thefunction. I found this post and tried making the suggested changes using do.call()

y=nit

X=X.mod

data <- list(y = nit, X=X.mod, n = n, l = l) 

parms <- c("beta","mu")

n.iter=3000

n.chains=2

n.burnin=1000

n.thin=5

samps <- do.call(jags.parallel,list(names(data),inits,parms,

                               "D:/Dropbox/PhD/Loading/Bayesian_Analysis/ExponentialGammaPriors.JAG",

                               n.iter,n.chains,n.burnin,n.thin))

# but then I get another error

Error in checkForRemoteErrors(lapply(cl, recvResult)) : 

  2 nodes produced errors; first error: Number of initialized chains (length(inits)) != n.chains

# I have 13 "beta" parameters and 721 "mu" parameters. I am confused why it seems to be asking for the number of chains to equal the number of initial estimates for all parameters? 

# if I set inits to something I get this error

Error in file(con, "w") : all connections are in use

Thank you for your help! 

[Carl Boettiger]

You need to initialize each chain (from ?jags.parallel: "a list with 'n.chains' elements").

You can specify a function that uses random starting values if you don't want to pick specific values for each chain.  I do recommend you see the examples in the documentation of ?jags.parallel, which include an illustration of such a function.  

p.s. good call with the `do.call`, that seems to be a common bug in the R2jags interface where it confuses variables for values if you don't use do.call...  And you got 2 nodes to produce errors simultaneously -- so that means something is happening in parallel!

Cheers,

[Katherine Lockhart]

Thanks Carl! I set up the initial values as a list, but I was still getting the connection error. I'm thinking R is opening connections and not closing them, so they are getting maxed out somehow in the process of connecting to multiple cores. 

I ultimately tried the code that Julie initially suggested using the wrapper function and the snow package to time the results and I can get it to run with no errors. The results are also nearly identical to when I run it on a single core. It's definitely faster for short MCMC chains and short burn in (2-5 times), however, the run times don't seem to be improving very much for very long model runs (like 10%). I'll keep testing it.  

Thanks again!