seed in RStan (vs Stan)

[Jiqiang Guo]

There are some difference between RStan and Stan in terms of how the seed is used. In RStan, I think taking a look at sampler_command in

https://code.google.com/p/stan/source/browse/rstan/rstan/inst/include/rstan/stan_fit.hpp

will help you.  I can remember there are at least tow differences:

1. initial values, in RStan, it is initialized once 

2. reset seed before sampling. (Line 475) 

So I think if you reset the seed in Stan before sampling, you might be able to get the same behavior. 

--

Jiqiang

On Wed, Oct 17, 2012 at 6:37 PM, Daniel Lee <bea...@alum.mit.edu> wrote:

On Wed, Oct 17, 2012 at 11:52 AM, Ben Goodrich <goodri...@gmail.com> wrote:
> On Wednesday, October 17, 2012 11:34:25 AM UTC-4, Daniel Lee wrote:
>>
>> Can you provide the model, data, and seed that causes the chain to go
>> off the rails? If there is a bug, perhaps we can hunt it down.
>
>
> https://groups.google.com/d/msg/stan-users/OjMSj1pITAQ/Tv8Pgy_88BoJ
> https://groups.google.com/d/msg/stan-dev/6hRAFUT1Pdk/Ye5OngvwQ74J
>
> Ben

Ben, can you help me figure out which seeds are causing these issues?
I see the problem by running rstan, I've pulled the seeds using
get_seed(chain[[3]]), but I'm not seeing the problem from the command
line.

--

[Daniel Lee]

On Wed, Oct 17, 2012 at 6:50 PM, Jiqiang Guo <guo...@gmail.com> wrote:
> There are some difference between RStan and Stan in terms of how the seed is
> used. In RStan, I think taking a look at sampler_command in
> https://code.google.com/p/stan/source/browse/rstan/rstan/inst/include/rstan/stan_fit.hpp
> will help you. I can remember there are at least tow differences:
> 1. initial values, in RStan, it is initialized once
> 2. reset seed before sampling. (Line 475)
>
> So I think if you reset the seed in Stan before sampling, you might be able
> to get the same behavior.
>

I didn't realize that the code was copied over. Bob added a few lines
of code that may solve a lot of those problems. I don't think we
realized that the change needed to be made in two places. See lines
346 to 369 in src/stan/gm/command.hpp.
https://code.google.com/p/stan/source/browse/src/stan/gm/command.hpp

Those lines should be included in stan_fit.hpp after 454.

If this code is being shared, I think we should pull out the common
code and maintain it in one place.



Daniel

[Jiqiang Guo]

I am going to include that and see if it solves that particular problem. 

--

Jiqiang

Daniel

--

[Ben Goodrich]

If you have a rstan rolled from master, you can do this in R

library(rstan)
example(stan_demo)
library(parallel)

models <- mclapply(1:83, mc.cores = 8, FUN = function(m) mean(unlist(get_logposterior(stan_demo(m, chains = 1, refresh = -1, seed = 12345), inc_warmup = TRUE)))) # better to use clang++ here

names(models) <- model_names

which(sapply(models, FUN = function(m) !is.na(m) && is.infinite(m)))

oxford  salm2   salm 

    32     36     37 

So, those three are suspicious. The oxford model throws many

Warning (non-fatal rejection): Error in function stan::prob::binomial_log(N4stan5agrad3varE): Probability parameter is -nan:0, but must be finite!

The salm and salm2 models throw many

Warning (non-fatal rejection): Error in function stan::prob::normal_log(N4stan5agrad3varE): Random variable is nan:0, but must not be nan!

In contrast, command-line Stan does not show a problem (each was run with --seed=12345).

goodrich@CYBERPOWERPC:/opt/stan/src/models/bugs_examples/vol1/oxford$ R --quiet -e "mean(read.csv('samples.csv', comment.char='#')\$lp__)"

> mean(read.csv('samples.csv', comment.char='#')$lp__)

[1] -4528.314

goodrich@CYBERPOWERPC:/opt/stan/src/models/bugs_examples/vol1/salm$ R --quiet -e "mean(read.csv('samples.csv', comment.char='#')\$lp__)"

> mean(read.csv('samples.csv', comment.char='#')$lp__)

[1] 1288.49

goodrich@CYBERPOWERPC:/opt/stan/src/models/bugs_examples/vol1/salm$ R --quiet -e "mean(read.csv('samples.csv', comment.char='#')\$lp__)"

> mean(read.csv('samples.csv', comment.char='#')$lp__)

[1] 1289.317

@Daniel, it would be nice if we could execute make SEED=12345 test-models to use that SEED always, which would make it easier to do these comparisons with rstan

@Mav. These models don't compile from rstan, so we should fix that somehow.

> which(sapply(models, is.na))

       rats_vec irt2_multilevel            irt2  irt_multilevel             irt 

             35              76              77              78              79 

       logistic 

             80 

Ben

[Ben Goodrich]

On Wednesday, October 17, 2012 11:39:29 PM UTC-4, Ben Goodrich wrote:

If you have a rstan rolled from master, you can do this in R

library(rstan)
example(stan_demo)
library(parallel)

models <- mclapply(1:83, mc.cores = 8, FUN = function(m) mean(unlist(get_logposterior(stan_demo(m, chains = 1, refresh = -1, seed = 12345), inc_warmup = TRUE)))) # better to use clang++ here

names(models) <- model_names

which(sapply(models, FUN = function(m) !is.na(m) && is.infinite(m)))

oxford  salm2   salm 

    32     36     37 

So, those three are suspicious.

That was before incorporating Mav's latest push. Now none of the average lp__ are infinite. And oxford and salm2 look reasonable. However, salm still has issues

> salm <- stan_demo("salm", seed = 12345)

TRANSLATING MODEL 'salm' FROM Stan CODE TO C++ CODE NOW.

COMPILING THE C++ CODE FOR MODEL 'salm' NOW.
[snip warnings]

1 warning generated.

SAMPLING FOR MODEL 'salm' NOW (CHAIN 1).

Iteration:    1 / 2000 [  0%]  (Adapting)

Iteration:  200 / 2000 [ 10%]  (Adapting)

Iteration:  400 / 2000 [ 20%]  (Adapting)

Iteration:  600 / 2000 [ 30%]  (Adapting)

Iteration:  800 / 2000 [ 40%]  (Adapting)

Iteration: 1000 / 2000 [ 50%]  (Adapting)

Iteration: 1200 / 2000 [ 60%]  (Sampling)

Iteration: 1400 / 2000 [ 70%]  (Sampling)

Iteration: 1600 / 2000 [ 80%]  (Sampling)

Iteration: 1800 / 2000 [ 90%]  (Sampling)

Iteration: 2000 / 2000 [100%]  (Sampling)

SAMPLING FOR MODEL 'salm' NOW (CHAIN 2).

Iteration:    1 / 2000 [  0%]  (Adapting)

Iteration:  200 / 2000 [ 10%]  (Adapting)

Iteration:  400 / 2000 [ 20%]  (Adapting)

Iteration:  600 / 2000 [ 30%]  (Adapting)

Iteration:  800 / 2000 [ 40%]  (Adapting)

Iteration: 1000 / 2000 [ 50%]  (Adapting)

Iteration: 1200 / 2000 [ 60%]  (Sampling)

Iteration: 1400 / 2000 [ 70%]  (Sampling)

Iteration: 1600 / 2000 [ 80%]  (Sampling)

Iteration: 1800 / 2000 [ 90%]  (Sampling)

Iteration: 2000 / 2000 [100%]  (Sampling)

SAMPLING FOR MODEL 'salm' NOW (CHAIN 3).

Iteration:    1 / 2000 [  0%]  (Adapting)

Iteration:  200 / 2000 [ 10%]  (Adapting)

Iteration:  400 / 2000 [ 20%]  (Adapting)

Iteration:  600 / 2000 [ 30%]  (Adapting)

Iteration:  800 / 2000 [ 40%]  (Adapting)

Iteration: 1000 / 2000 [ 50%]  (Adapting)

Iteration: 1200 / 2000 [ 60%]  (Sampling)

Iteration: 1400 / 2000 [ 70%]  (Sampling)

Iteration: 1600 / 2000 [ 80%]  (Sampling)

Iteration: 1800 / 2000 [ 90%]  (Sampling)

Iteration: 2000 / 2000 [100%]  (Sampling)

SAMPLING FOR MODEL 'salm' NOW (CHAIN 4).

Error in sampler$call_sampler(args_list[[i]]) : 

  No acceptably small step size could be found. Perhaps the posterior is not continuous?

error occurred during calling the sampler; sampling not done

Ben

[Daniel Lee]

> @Daniel, it would be nice if we could execute make SEED=12345 test-models to
> use that SEED always, which would make it easier to do these comparisons
> with rstan

Sorry, that's pretty hard to do. The models aren't run from the
makefile. They are run from within google test using popen().

The good news: if you've already run test-models, you can just run the
model using the seed you want. The executables are all built and ready
to go.

[Bob Carpenter]

On 10/17/12 7:36 PM, Daniel Lee wrote:
> On Wed, Oct 17, 2012 at 6:50 PM, Jiqiang Guo <guo...@gmail.com> wrote:
>> There are some difference between RStan and Stan in terms of how the seed is

>> ...

> I didn't realize that the code was copied over.

I keep forgetting, too, though I know it has to
be done that way (separate commands).

> Bob added a few lines
> of code that may solve a lot of those problems. I don't think we
> realized that the change needed to be made in two places. See lines
> 346 to 369 in src/stan/gm/command.hpp.
> https://code.google.com/p/stan/source/browse/src/stan/gm/command.hpp
>
> Those lines should be included in stan_fit.hpp after 454.
>
> If this code is being shared, I think we should pull out the common
> code and maintain it in one place.

That's a great idea. Maybe the best practical course is
to have the three of us (me, Daniel, Jiqiang) sit down and
review what's shared and how we might consolidate it.

- Bob

>
> Daniel
>