Parallel chains now supported (sort of)

[Ben Goodrich]

Good news everyone (particularly people who are willing to help us test)! Jiqiang wrote a function that will be included in rstan 1.0.3 that facilitates running chains in parallel. In particular, you can do something like

stopifnot(require(rstan))

schools_dat <- list(J = 8,
                    y = c(28,  8, -3,  7, -1,  1, 18, 12),
                    sigma = c(15, 10, 16, 11,  9, 11, 10, 18))

eight_schools <- stan(file = "http://wiki.stan.googlecode.com/git/rstangettingstarted/8schools.stan",
                      data = schools_dat, iter = 1, chains = 1) # just get it to compile the model

library(parallel) # or some other parallelizing package
chains <- 8

posterior <- mclapply(1:chains, mc.cores = chains, FUN = function(chain) { # or another parallelizing function
               stan(fit = eight_schools, data = schools_dat, chains = 1, verbose = FALSE, refresh = -1, chain_id=chain)
             })

posterior <- sflist2stanfit(posterior) # see below for function definition
print(posterior)
plot(posterior)

where the necessary new function sflist2stanfit() is defined at the end of this post and is also attached. 

Some things to note. In the example above, I first ran the eight schools model for one iteration of one chain just to get the code to compile. Then I passed that eight_schools object to the fit argument of stan() in each of the eight chains (with the defaults for warmup, iter, etc.). This route avoids having to recompile the code eight times and will work for mclapply() in the parallel package, which executes the function locally on a multicore machine. Note that mclapply() does not do anything in parallel on Windows, so you need to use an alternative parallelizing library / function in that case.

This route will probably not work if you are using the snow package or a similar package to execute the chains on a remote server because the remote server will not know what to do with the code compiled on your laptop. In that case, you would need to specify file or model_code rather than fit, so that the remote server would compile the code for each chain. Also, if file is a filepath on your laptop, the remote server will not be able to access it, so you would need to upload your .stan file to a readable directory on the remote server and specify file accordingly. Similarly, the data somehow need to be available on the remote server and the remote server obviously needs to have rstan installed and loaded. Note that stanfit objects can be quite large in size, so you could be pulling GB of uncompressed data from the remote server to your laptop at the end.

If you have access to some sort of computer lab with a bunch of physical PCs that have rstan but are not networked, you could execute stan() with chains = 1 on each PC, save the  resulting stanfit object to a USB stick with something like save(chain_3, file = "chain_3.RData", compress = "xz"), use the load() command to load each .RData object onto a common machine, put them into a list, and use sflist2stanfit() to consolidate them into a stanfit object that can be analyzed by the tools in the rstan package. Note that each chain needs to be given a different object name to avoid them overwriting each other when you load() them onto the common machine.

For the high-performance (including parallel) capabilities that are available in R or R packages, execute help(package = parallel) in R or see 

http://cran.r-project.org/web/views/HighPerformanceComputing.html

Ben (but mostly Jiqiang)

P.S. The new function below can be pasted in now but is now defined in the rstan tree at the end of

https://code.google.com/p/stan/source/browse/rstan/rstan/R/stanfit-class.R

sflist2stanfit <- function(sflist) {
  # merge a list of stanfit objects into one
  # Args:
  #   sflist, a list of stanfit objects, each element of which 
  #   should have equal length of `iter`, `warmup`, and `thin`. 
  # Returns: 
  #   A new stanfit objects have all the chains in each element of sf_list. 
  #   The date would be where the new object is created. 
  # Note: 
  #   * method get_seed would not work well for this merged 
  #     stanfit object. But all the information is still there. 
  #   * When print function is called, the sampler is obtained 
  #     only from the first chain. 
  #  
  sf_len <- length(sflist) 
  if (sf_len < 2) stop("'sflist' should have more than 1 elements")
  if (!is.list(sflist) || 
      any(sapply(sflist, function(x) !is(x, "stanfit")))) {   
    stop("'sflist' must be a list of 'stanfit' objects")
  }
  if (any(sapply(sflist, function(x) x@mode != 0))) {
    stop("each 'stanfit' object in 'sflist' must contain samples") 
  }   
  for (i in 2:sf_len) { 
    if (!identical(sflist[[i]]@sim$pars_oi, sflist[[1]]@sim$pars_oi) || 
        !identical(sflist[[i]]@sim$dims_oi, sflist[[1]]@sim$dims_oi))   
      stop("parameters in element ", i, " (stanfit object) is different from in element 1") 
    if (sflist[[i]]@sim$n_save[1] != sflist[[1]]@sim$n_save[1] ||  
        sflist[[i]]@sim$warmup2[1] != sflist[[1]]@sim$warmup2[1])
      stop("all 'stanfit' objects should have equal length of iterations and warmup") 
  } 
  n_chains = sum(sapply(sflist, function(x) x@sim$chains))  
  sim = list(samples = do.call(c, lapply(sflist, function(x) x@sim$samples)),  
             chains = n_chains, 
             iter = sflist[[1]]@sim$iter, 
             thin = sflist[[1]]@sim$thin, 
             warmup = sflist[[1]]@sim$warmup,
             n_save = rep(sflist[[1]]@sim$n_save[1], n_chains), 
             warmup2 = rep(sflist[[1]]@sim$warmup2[1], n_chains), 
             thin = rep(sflist[[1]]@sim$thin[1], n_chains), 
             permutation = do.call(c, lapply(sflist, function(x) x@sim$permutation)), 
             pars_oi = sflist[[1]]@sim$pars_oi, 
             dims_oi = sflist[[1]]@sim$dims_oi, 
             fnames_oi = sflist[[1]]@sim$fnames_oi,
             n_flatnames = sflist[[1]]@sim$n_flatnames) 
  nfit <- new("stanfit", 
              model_name = sflist[[1]]@model_name, 
              model_pars = sflist[[1]]@model_pars,
              par_dims  = sflist[[1]]@par_dims, 
              mode = 0L, 
              sim = sim, 
              inits = do.call(c, lapply(sflist, function(x) x@inits)), 
              stan_args = do.call(c, lapply(sflist, function(x) x@stan_args)), 
              stanmodel = sflist[[1]]@stanmodel, 
              date = date(), 
              .MISC = new.env()) 
  invisible(nfit)
} 

[Andrew Gelman]

> install.packages("parallel")

Installing package(s) into ‘/Library/Frameworks/R.framework/Versions/2.15/Resources/library’

(as ‘lib’ is unspecified)

Warning in install.packages :

  package ‘parallel’ is not available (for R version 2.15.1)

On Oct 12, 2012, at 12:14 AM, Ben Goodrich wrote:

Good news everyone (particularly people who are willing to help us test)! Jiqiang wrote a function that will be included in rstan 1.0.3 that facilitates running chains in parallel. In particular, you can do something like

stopifnot(require(rstan))

schools_dat <- list(J = 8,
                    y = c(28,  8, -3,  7, -1,  1, 18, 12),
                    sigma = c(15, 10, 16, 11,  9, 11, 10, 18))

eight_schools <- stan(file = "http://wiki.stan.googlecode.com/git/rstangettingstarted/8schools.stan",
                      data = schools_dat, iter = 1, chains = 1) # just get it to compile the model

library(

parallel) # or some other parallelizing package

chains <- 8

posterior <- mclapply(1:chains, mc.cores = chains, FUN = function(chain) { # or another parallelizing function

               stan(fit = eight_schools, data = schools_dat, chains = 1, verbose = FALSE, refresh = -1)
             })

posterior <- sflist2stanfit(posterior) # see below for function definition
print(posterior)
plot(posterior)

--
 
 
<sflist2stanfit.txt>

[Ben Goodrich]

On Friday, October 12, 2012 12:18:12 AM UTC-4, Andrew Gelman wrote:

> install.packages("parallel")

Installing package(s) into ‘/Library/Frameworks/R.framework/Versions/2.15/Resources/library’

(as ‘lib’ is unspecified)

Warning in install.packages :

  package ‘parallel’ is not available (for R version 2.15.1)

The parallel package comes with recent versions of R, so you shouldn't need to install it separately. Just

library(parallel)

in R.

Ben

[Andrew Gelman]

Yup.

--
 
 

[Sebastian Weber]

Hi!

For the lazy people out there, here is a convenient function that wraps up things from Andrew:

stanParallel <- function(fit=NULL, chains=4, ...) {
args <- list(...)
args$chains <- 1

## if fit is given, then use the compiled code from there
if(is.null(fit)) {
targs <- args
targs$iter <- 1
fit <- do.call(stan, targs)
}

args$file <- NULL
args$fit <- fit
args$refresh <- -1
args$verbose <- FALSE

posterior <- mclapply(1:chains, FUN = function(chain) do.call(stan, args) )

posterior <- sflist2stanfit(posterior) # see below for function definition

invisible(posterior)
}

Just worked like a charm for me.

One side note: You have to use the mc.cores option to set the number of cores to use, i.e. put a options(mc.cores=2) in your .Rprofile .

Best,

Sebastian

> --
>
>

--
Sebastian Weber
Group of Cell Communication and Control
Freiburg Institute for Advanced Studies - FRIAS
School of Life Sciences - LIFENET
Albert-Ludwigs-Universität Freiburg
Albertstr. 19
79104 Freiburg
T.: +49-761-203-97237
Fax:+49-761-203-97334

[Ben Goodrich]

I failed to mention something potentially important. In code such this:

On Friday, October 12, 2012 2:41:30 AM UTC-4, Sebastian Weber wrote:

    posterior <- mclapply(1:chains, FUN = function(chain) do.call(stan, args) )

It is important to pass chain_id=chain to Stan, as in

posterior <- mclapply(1:chains, FUN = function(chain) {
               args$chain_id <- chain
               do.call(stan, args)
             })

 

rstan does not use R's pseduo random number generators; rather Stan's uses one of Boost's PRNGs. I'm still not convinced that what we are doing is optimal, but setting chain_id causes Stan to make some effort to ensure that the chains start at interspersed points in the period of the PRNG.

Ben
 

[Sebastian Weber]

Thanks for noting, but your code has a potential concurrency problem, if I am right.

I would opt for this line of code

posterior <- mclapply(1:chains, FUN = function(chain){ do.call(stan, c(args, list(chain_id=chain))) } )

This avoids that args$chain_id gets overwritten by using a local copy in any case, but this is maybe overly cautious…

Best,

Sebastian

[Bob Carpenter]

On 10/12/12 4:28 AM, Ben Goodrich wrote:
> I failed to mention something potentially important. In code such this:
>
> On Friday, October 12, 2012 2:41:30 AM UTC-4, Sebastian Weber wrote:
>
> posterior <- mclapply(1:chains, FUN = function(chain) do.call(stan, args) )
>
>
> It is important to pass chain_id=chain to Stan, as in
>
> |
> posterior <- mclapply(1:chains, FUN = function(chain) {
> args$chain_id <- chain
> do.call(stan, args)
> })
> |
>
> rstan does not use R's pseduo random number generators; rather Stan's uses one of Boost's PRNGs. I'm still not convinced

> that what we are doing is optimal,but setting chain_id causes Stan to make some effort to ensure that the chains start

> at interspersed points in the period of the PRNG.

Optimal in what sense? It's very easy to change.

I should explain what we're doing more clearly in
the manual. I plugged in the L'Ecuyer RNG with
a seed for chain 1. Then each subsequent chain
uses skip-ahead so the random numbers don't clash.
I used a very very large skip ahead that would
require a whole lot of random numbers per chain.
I'll look it up and put the right numbers in.

The period of L'Ecuyer is pretty good, but
speed isn't optimal. And if you do things like
use consecutive or close seeds, you get closely
related random numbers at least initially.

The C++ code lets you plug and play Boost RNGs
(or other RNGs implementing their concept.)

- Bob

>
> Ben
>
> --
>
>

[Jiqiang Guo]

I think what Ben was trying to say it is *not* possible to pass chain_id to stan().  Whatever you give as chain_id for stan() (as it has ... argument, so it would not complain), the real chain ids would be from 1 to the number of chains (parameter chains).  So it is important to have very different seeds for stan if run in parallel in rstan 1.0.2. 

So I should let stan() honor user supplied chain_id in rstan 1.0.3.  Then the best way would be to let all stan call have the same seed, but different chain_id's. 

--

Jiqiang

--

[Tim Triche, Jr.]

you are a hero Jiqiang, thanks for your many efforts on supporting stan and rStan

--
 
 

--
A model is a lie that helps you see the truth.

Howard Skipper

[Ben Goodrich]

On Friday, October 12, 2012 9:49:03 AM UTC-4, Jiqiang Guo wrote:

I think what Ben was trying to say it is *not* possible to pass chain_id to stan().  Whatever you give as chain_id for stan() (as it has ... argument, so it would not complain), the real chain ids would be from 1 to the number of chains (parameter chains).  

Ah. I thought stan() would pass everything in its ... to Stan, but maybe not until 1.0.3. In any event, supplying the chain_id does not hurt now, and it may do something good in a future version of rstan / Stan.
 

So it is important to have very different seeds for stan if run in parallel in rstan 1.0.2. 

As it stands right now, using mclapply() is not so bad because stan() by default passes a pseudo-random seed from R to Stan because

seed = sample.int(.Machine$integer.max, 1)

by default. Those $chains values are realized sequentially in the main R process before it gets forked $chains times and starts $chains Stan processes. That is not the best thing in the world, but if you are only doing a few chains for a few thousand iterations then it probably won't compromise the independence of the chains too badly.

If you were using another parallelizing function (such as on Windows) that does not fork the main R process but actually starts $chains new R processes in the background at approximately the same time, then sample.int() could produce the same or similar seeds. I was hoping that specifying chain_id would help, but I guess rstan 1.0.2 doesn't pass that to Stan.
 

So I should let stan() honor user supplied chain_id in rstan 1.0.3.  Then the best way would be to let all stan call have the same seed, but different chain_id's. 

That would be good. Particularly if we upgrade from the Boost PRNG we are using now.

Ben

[atzap]

For us Windows folks, it doesn't look like the parallel package works.  Here is my attempt using snow.

# compile model

fit.stan <- stan(model_code = m2.SEM, data = bdat, iter = 1, chains = 1, verbose = T)

# load snow

library(snow)

# setup cluster with two cores

cl <- makeCluster(c("localhost","localhost"), type = "SOCK")

# create list with data, parameters to be traced, compiled model, and function for generating initial values

cldata <- list(data=bdat, params=params, fit = fit.stan, inits=its)

clusterEvalQ(cl, library(rstan))    # loads rjags library on each node

clusterExport(cl, "cldata")         # loads cldata on each node

stanmod <- function(i) {

                        stan(fit = cldata$fit, data = cldata$data, inits = cldata$inits(),

                                  iter = 700, warmup = 250, thin = 2, chains = 1,

                                  pars = cldata$params)

}

stan.out <- clusterApply(cl, 1:2, stanmod)

stopCluster(cl)

Obviously, this function can be cleaned up to supply the iter, warmup, and thin values directly to stanmod.  

More importantly, you probably want to follow the sflist2stanfit approach and pull out the samples for post-processing using stan.out[[1]]@sim$samples or something like as.mcmc.list(lapply(stan.out, function(x){as.mcmc(extract(x, inc_warmup = F)[,1,])})).

Hope this isn't too daft!

Andrew

[Ben Goodrich]

On Friday, October 12, 2012 2:41:09 PM UTC-4, atzap wrote:

For us Windows folks, it doesn't look like the parallel package works.  Here is my attempt using snow.

stanmod <- function(i) {

                        stan(fit = cldata$fit, data = cldata$data, inits = cldata$inits(),

                                  iter = 700, warmup = 250, thin = 2, chains = 1,

                                  pars = cldata$params)

}

stan.out <- clusterApply(cl, 1:2, stanmod)

stopCluster(cl)

Here, you would want to be careful about seeds because in your example snow will start two R processes locally at essentially the same time. Apparently, recent versions of R initiate the PRNG state with some function of the time and the process ID. Maybe the process ID sufficiently distinguishes the states. But if you have the random package installed and a live internet connection, then you could add

seed = random:::randomNumbers(1, min = -10^9, max = 10^9, col = 1)[1,1] + 10^9 + 1

to the stan() call to download one random integer per chain from www.random.org (which uses atmospheric noise to create randomness).

Ben

[Ben Goodrich]

On Friday, October 12, 2012 4:28:25 AM UTC-4, Ben Goodrich wrote:

I failed to mention something potentially important. In code such this:

On Friday, October 12, 2012 2:41:30 AM UTC-4, Sebastian Weber wrote:

    posterior <- mclapply(1:chains, FUN = function(chain) do.call(stan, args) )

It is important to pass chain_id=chain to Stan, as in

You can also use the chain argument to switch whether to print the progress to the screen

posterior <- mclapply(1:chains, FUN = function(chain) { # or another parallelizing function

  stan(fit = eight_schools, data = schools_dat, chains = 1, verbose = FALSE,

       refresh = ifelse(chain == 1, 200, -1), chain_id = chain)
})

However, this only works properly when R is running in a shell (not RStudio or various other GUIs that hold the output back until all chains are done). And the other chains may be progressing at a substantially faster or slower rate than the first chain.

Ben

[W Robert Long]

Hi all

Is anyone running Stan on CentOS linux ? Do you anticipate any issues ?
I should have a new quad-core PC next week with CentOS 5.4 as a dual
boot and hope to play around with the parallel chains.

Thanks
Rob

[Marcus Brubaker]

I routinely run Stan on CentOS 6 but I have never used R or RStan.  Before I upgraded I had done a little work with 5.4 but I'm not sure how it would run currently.  The only issue I would expect to come up is that g++ is old and will be slow to compile on 5.4.  Depending just how old it is, it may also choke a little on Stan's heavy use of templates through Boost and Eigen.

Cheers,
Marcus

--

[Ben Goodrich]

In general, there is no reason why Stan should not work on any Linux distribution that has the dependencies, which are all widely-used and should be available in binary form from the Linux distribution's package repositories.

Ben

[Song]

I got the following error when running the 8 school examples on a MacBook:

posterior <- sflist2stanfit(posterior) # see below for function definition

Error in FUN(X[[1L]], ...) : 

  no slot of name "mode" for this object of class "stanfit"

Song

On Friday, October 12, 2012 12:14:45 AM UTC-4, Ben Goodrich wrote:

Good news everyone (particularly people who are willing to help us test)! Jiqiang wrote a function that will be included in rstan 1.0.3 that facilitates running chains in parallel. In particular, you can do something like

stopifnot(require(rstan))

schools_dat <- list(J = 8,
                    y = c(28,  8, -3,  7, -1,  1, 18, 12),
                    sigma = c(15, 10, 16, 11,  9, 11, 10, 18))

eight_schools <- stan(file = "http://wiki.stan.googlecode.com/git/rstangettingstarted/8schools.stan",
                      data = schools_dat, iter = 1, chains = 1) # just get it to compile the model

library(

parallel) # or some other parallelizing package
chains <- 8

posterior <- mclapply(1:chains, mc.cores = chains, FUN = function(chain) { # or another parallelizing function
               stan(fit = eight_schools, data = schools_dat, chains = 1, verbose = FALSE, refresh = -1, chain_id=chain)
             })

posterior <- sflist2stanfit(posterior) # see below for function definition
print(posterior)
plot(posterior)

[Jiqiang Guo]

Please see comments below. 

On Thu, Oct 18, 2012 at 11:45 AM, Song <mdq...@gmail.com> wrote:

I got the following error when running the 8 school examples on a MacBook:

posterior <- sflist2stanfit(posterior) # see below for function definition

Error in FUN(X[[1L]], ...) : 

  no slot of name "mode" for this object of class "stanfit"

The reason might be that you are not using rstan 1.0.2.  Please report at least the version 

number of rstan when reporting problems.  Thanks. 

 

--
 
 

[Song]

You are right, 1.0.1 was still on that machine.

[Song]

You may also want a line like this:

  targs$thni <- 1

to avoid error when thin is supplied in the ...

On Friday, October 12, 2012 2:41:30 AM UTC-4, Sebastian Weber wrote:

[Tom Wallis]

Hi everyone,

I'm resurrecting this thread to ask whether people have an opinion on the use of an mc.preschedule=FALSE flag in the mclapply function.

From the help there:

if set to TRUE then the computation is first divided to (at most) as many jobs are there are cores and then the jobs are started, each job possibly covering more than one value. If set to FALSE then one job is forked for each value of X. The former is better for short computations or large number of values in X, the latter is better for jobs that have high variance of completion time and not too many values of X compared to mc.cores.

It seems to me that the latter is going to be the situation corresponding best to Stan. Has anyone compared these strategies?

Thanks

Tom

On Friday, October 12, 2012 6:14:45 AM UTC+2, Ben Goodrich wrote:

Good news everyone (particularly people who are willing to help us test)! Jiqiang wrote a function that will be included in rstan 1.0.3 that facilitates running chains in parallel. In particular, you can do something like

stopifnot(require(rstan))

schools_dat <- list(J = 8,
                    y = c(28,  8, -3,  7, -1,  1, 18, 12),
                    sigma = c(15, 10, 16, 11,  9, 11, 10, 18))

eight_schools <- stan(file = "http://wiki.stan.googlecode.com/git/rstangettingstarted/8schools.stan",
                      data = schools_dat, iter = 1, chains = 1) # just get it to compile the model

library(parallel) # or some other parallelizing package

chains <- 8

posterior <- mclapply(1:chains, mc.cores = chains, FUN = function(chain) { # or another parallelizing function
               stan(fit = eight_schools, data = schools_dat, chains = 1, verbose = FALSE, refresh = -1, chain_id=chain)

             })

posterior <- sflist2stanfit(posterior) # see below for function definition
print(posterior)
plot(posterior)

[Bob Carpenter]

Sounds like an empirical question that may also vary by hardware.
So the best thing to do would be to try it on whatever setup you're
going to use. We'd love to hear the results.

You probably don't want to fork more processes than you have cores,
because then the OS has to do the swapping of processes over cores.
This probably won't be too bad unless you wind up also having to swap
large chunks of memory.

I'm not an R programmer, but the standard solution to this is some
kind of execution queue. You push all the processes that need to be run
onto a queue. Then, for each core you want to use, you pull a process
off the queue. If a core's process finishes, it pulls the next job.
Scheduling processes across core is, in general, a complicated business,
especially when the process resource usage and execution priority varies.

Hyperthreading on Intel machines makes tuning confusing. A CPU with
hyperthreading reports twice as many cores as there are physical cores.
My own experience is to cap the number of processes at the actual number
of cores, not the number of virtual cores, but your mileage may vary.

- Bob

> to get the code to compile. Then I passed that eight_schools object to the _fit_ argument of _stan()_ in each of the

> eight chains (with the defaults for warmup, iter, etc.). This route avoids having to recompile the code eight times

> and will work for _mclapply()_ in the parallel package, which executes the function locally on a multicore machine.
> Note that _mclapply()_ does not do anything in parallel on Windows, so you need to use an alternative parallelizing

> library / function in that case.
>
> This route will probably not work if you are using the snow package or a similar package to execute the chains on a

> *remote server* because the remote server will not know what to do with the code compiled on your laptop. In that
> case, you would need to specify _file_ or _model_code_ rather than _fit_, so that the remote server would compile
> the code for each chain. Also, if _file_ is a filepath on your laptop, the remote server will not be able to access

> it, so you would need to upload your .stan file to a readable directory on the remote server and specify

> _file_ accordingly. Similarly, the data somehow need to be available on the remote server and the remote server

> obviously needs to have rstan installed and loaded. Note that stanfit objects can be quite large in size, so you
> could be pulling GB of uncompressed data from the remote server to your laptop at the end.
>
> If you have access to some sort of computer lab with a bunch of physical PCs that have rstan but are not networked,

> you could execute _stan()_ with _chains = 1_ on each PC, save the resulting stanfit object to a USB stick with
> something like _save(chain_3, file = "chain_3.RData", compress = "xz")_, use the _load()_ command to load each
> .RData object onto a common machine, put them into a list, and use _sflist2stanfit()_ to consolidate them into a

> stanfit object that can be analyzed by the tools in the rstan package. Note that each chain needs to be given a

> different object name to avoid them overwriting each other when you _load()_ them onto the common machine.

>
> For the high-performance (including parallel) capabilities that are available in R or R packages,

> _execute help(package = parallel)_ in R or see

> --
> You received this message because you are subscribed to the Google Groups "stan users mailing list" group.
> To unsubscribe from this group and stop receiving emails from it, send an email to stan-users+...@googlegroups.com.
> For more options, visit https://groups.google.com/groups/opt_out.

[Katherine Grace Jonas]

On Oct 12, 2012, at 12:14 AM, Ben Goodrich wrote:

Good news everyone (particularly people who are willing to help us test)! Jiqiang wrote a function that will be included in rstan 1.0.3 that facilitates running chains in parallel. In particular, you can do something like

stopifnot(require(rstan))

schools_dat <- list(J = 8,
                    y = c(28,  8, -3,  7, -1,  1, 18, 12),
                    sigma = c(15, 10, 16, 11,  9, 11, 10, 18))

eight_schools <- stan(file = "http://wiki.stan.googlecode.com/git/rstangettingstarted/8schools.stan",
                      data = schools_dat, iter = 1, chains = 1) # just get it to compile the model

When using this method to sample in parallel, does each chain use the same initial values?  If the initial values are set in running this single iteration, and that iteration is passed to the parallel operation, would each chain would follow from the same set of initial values?  What if, in the example above, iter=0?

[Ben Goodrich]

On Thursday, July 3, 2014 5:27:10 PM UTC-4, Katherine Grace Jonas wrote:

When using this method to sample in parallel, does each chain use the same initial values?  If the initial values are set in running this single iteration, and that iteration is passed to the parallel operation, would each chain would follow from the same set of initial values?  What if, in the example above, iter=0?

Regardless of whether you run the chains in parallel or serially, (R)Stan does not "continue" a previous run, even if you pass a stanfit object to the fit argument of stan(). In other words, each chain starts fresh currently. [I think we are working on having the ability to continue a previous run, but it is not much of a priority].

By default, (R)Stan randomly draws initial values. Even if you specify a seed (which you should probably do), as long as you specify the chain_id (which you should definitely do), then the state of the pseudo-random number generator will be different for each chain because it is advanced by a function of the chain_id.

If you wanted to specify initial values for each chain (rarely is this necessary), then you could do something like

initial_values <- list(something, different, for, each, chain)
posterior <- sflist2stanfit(mclapply(1:5, mc.cores = 5, FUN = function(i) {
  stan(fit = eight_schools, seed = 12345, chain_id = i, chains = 1, inits = initial_values[[i]])
}))

Ben

 

[Katherine Grace Jonas]

So the fit argument is taking the compiled model, but not the samples generated from it.  That makes sense.  Thank you Ben.  

[Bob Carpenter]

The stanfit object has both the samples and the
compiled model.

The only reason we don't support restarting is that
the C++ interface to Stan doesn't provide a way to set
the mass matrix and step size for NUTS. It also
doesn't provide a good way to continue with the same PRNG
state.

We do plan to enhance the C++ interface to allow proper
restarting (though we still haven't decided what to do about
restarting the PRNG).

This is among the things we're going to be updating,
but as Ben says, it's not a huge priority, so we may
not get to it very soon.

- Bob

> --
> You received this message because you are subscribed to the Google Groups "Stan users mailing list" group.

> To unsubscribe from this group and stop receiving emails from it, send an email to stan-users+...@googlegroups.com.

> For more options, visit https://groups.google.com/d/optout.

[Katherine Grace Jonas]

That's interesting.  Thanks for explaining.

[Marcus Brubaker]

As a side note (and perhaps one more for stan-dev), there are some easy ways to save/restart the PRNG.  If nothing else, one way is to reseed the PRNG before every iteration using a randomly generated seed.  This seed is generated at the end of the last iteration using the PRNG and can be saved with the state, making it easy to restart things.

Two potential downsides of this is it 1) could be a little slow depending on whether seeding is expensive for a particular PRNG and 2) you are potentially not using the entire period of the PRNG depending on how it works.  Of course for 1) it's unlikely to be so slow that we'd notice it for anything but the fastest models and for 2) I'm not sure this is really an issue but I could be wrong.

On the up side, it works for any PRNG, is simple to implement and is compatible with user-specified seeds.

Cheers,

Marcus