Fwd: [stan-users] Re: t-test (effects of initialization)

[Ben Goodrich]

It seems in the original example, the problem was basically a small sample size. But even if I multiple the sample sizes by 10 and add weakly informative priors (attached), the sampler goes off the rails

> print(fit, digits = 3)
Inference for Stan model: ch7_2.
4 chains: each with iter=21000; warmup=1000; thin=1; 21000 iterations saved.

            mean se_mean      sd      2.5%       25%       50%       75%     97.5% n_eff    Rhat
alpha1   104.930   2.214  17.274   104.533   105.139   105.224   105.309   151.083    61   1.001
alpha2    57.667  24.276  34.332    -1.946    57.299    77.431    77.546    77.738     2 308.550
sigma1     6.446   4.400  16.274     2.619     2.722     2.781     2.839    51.303    14   1.279
sigma2    21.894  23.489  33.242     2.548     2.676     2.758    19.808    83.114     2  25.105
dif      -15.449  25.262  38.869   -80.330   -17.096     0.014     0.137    48.565     2   3.143
lp__   -1978.780 375.129 668.558 -3276.437 -2866.694 -1521.770 -1520.283 -1519.298     3   2.412

On Tuesday, October 16, 2012 3:47:02 PM UTC-4, Bob Carpenter wrote:

On 10/16/12 3:32 PM, Shige Song wrote:
> Thanks for the help, Ben. Dropping the seed did lead to closer estimates of the sigma.

There's always a seed.  If you don't specify one, Stan
generates it randomly given the time.

> I did not expect setting random
> number seed can have such an effect.

It comes down to the initialization.  If you do
your own initial points, it won't matter nearly
as much for sampling.  But random inits, especially
with weak constraints, can lead to very bad initial
points that take a lot of time to adapt out of.

If you run long enough, adaptation should find a
better place and adapt to that before sampling.

- Bob

[Ben Goodrich]

This one still looks problematic for rstan as of commit 0cde41e but still looks okay from command-line Stan.

Ben

[Jiqiang Guo]

But I did not have problem with it in rstan as of commit a6ca1191aeac. 

Actually, if I comment out line 500 of 

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

RStan should have exactly the same result as Stan for the same seed. 

--

Jiqiang

> library(rstan)

Loading required package: Rcpp

Loading required package: inline

Loading required package: RcppEigen

Loading required package: Matrix

Loading required package: lattice

rstan (Version 1.0.2, packaged: 2012-10-18 14:46:14 UTC, GitRev: a6ca1191aeac)

> 

> set.seed(12345)

> 

> n1 <- 600                # Number of females

> n2 <- 400                # Number of males

> mu1 <- 105                # Population mean of females

> mu2 <- 77.5                # Population mean of males

> sigma <- 2.75                # Average population SD of both

> 

> n <- n1+n2                # Total sample size

> y1 <- rnorm(n1, mu1, sigma)        # Data for females

> y2 <- rnorm(n2, mu2, sigma)        # Date for males

> y <- c(y1, y2)                # Aggregate both data sets

> x <- rep(c(0,1), c(n1, n2))        # Indicator for male

> 

> ch7_2_data <- list(n1=n1, n2=n2, x=x, y1=y1, y2=y2)

> 

> fit <- stan(file="ch7_2.stan", data=ch7_2_data, warmup=1000, iter=21000)

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

COMPILING THE C++ CODE FOR MODEL 'ch7_2' NOW.

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

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

......

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

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

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

......

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

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

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

......

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

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

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

......

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

> 

> print(fit, digits=3)

Inference for Stan model: ch7_2.

4 chains: each with iter=21000; warmup=1000; thin=1; 21000 iterations saved.

            mean se_mean      sd      2.5%       25%       50%       75%

alpha1   104.856   0.458   4.297   104.951   105.142   105.218   105.296

alpha2    77.487   0.004   0.137    77.222    77.396    77.488    77.578

sigma1     3.158   0.485   4.273     2.619     2.715     2.767     2.824

sigma2     2.720   0.003   0.094     2.546     2.655     2.718     2.782

dif        0.437   0.487   4.276    -0.192    -0.034     0.049     0.133

lp__   -1537.059  20.151 147.852 -1526.114 -1521.783 -1520.692 -1519.952

           97.5% n_eff  Rhat

alpha1   105.443    88 1.030

alpha2    77.758  1209 1.009

sigma1     2.976    78 1.034

sigma2     2.912   775 1.002

dif        0.347    77 1.034

lp__   -1519.239    54 1.050

Samples were drawn using NUTS2 at Thu Oct 18 11:25:40 2012.

For each parameter, n_eff is a crude measure of effective sample size,

and Rhat is the potential scale reduction factor on split chains (at 

convergence, Rhat=1).

> 

--
 
 

[Ben Goodrich]

On Thursday, October 18, 2012 11:30:20 AM UTC-4, Jiqiang Guo wrote:

But I did not have problem with it in rstan as of commit a6ca1191aeac. 

Actually, if I comment out line 500 of 

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

RStan should have exactly the same result as Stan for the same seed. 

I think RStan should give exactly the same result as Stan for the same seed. But to me, the output still looks wrong. We shouldn't have a n_eff of less than 80 when 20,000 iterations are kept in a Gaussian model with unknown mean and variance and sample sizes of 600 and 400. When I do it from command-line Stan and look at the Naive SE of sigma1 in coda, it is only .0026.

Ben

[Bob Carpenter]

On 10/18/12 12:17 PM, Ben Goodrich wrote:
> On Thursday, October 18, 2012 11:30:20 AM UTC-4, Jiqiang Guo wrote:
>
> But I did not have problem with it in rstan as of commit a6ca1191aeac.
> Actually, if I comment out line 500 of
> https://code.google.com/p/stan/source/browse/rstan/rstan/inst/include/rstan/stan_fit.hpp

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

> RStan should have exactly the same result as Stan for the same seed.
>
>
> I think RStan should give exactly the same result as Stan for the same seed.

How could it not? Is the seed somehow not getting passed?

Is it because of something like different optimization
levels somehow mucking things up?

> But to me, the output still looks wrong. We
> shouldn't have a n_eff of less than 80 when 20,000 iterations are kept in a Gaussian model with unknown mean and
> variance and sample sizes of 600 and 400.

That does seem awfully low.

Does it have something to do with the difference
in the way n_eff is calculated? Coda uses some
spectral thing and RStan uses the split R-hat var+
estimator.

Can you see if the samples themselves are the same?

> When I do it from command-line Stan and look at the Naive SE of sigma1 in
> coda, it is only .0026.

Call me "Goldilocks", but that almost seems too good.

- Bob

[Ben Goodrich]

On Thursday, October 18, 2012 1:18:43 PM UTC-4, Bob Carpenter wrote:

On 10/18/12 12:17 PM, Ben Goodrich wrote:
> On Thursday, October 18, 2012 11:30:20 AM UTC-4, Jiqiang Guo wrote:
>
>     But I did not have problem with it in rstan as of commit a6ca1191aeac.
>     Actually, if I comment out line 500 of
>     https://code.google.com/p/stan/source/browse/rstan/rstan/inst/include/rstan/stan_fit.hpp
>     <https://code.google.com/p/stan/source/browse/rstan/rstan/inst/include/rstan/stan_fit.hpp>,
>     RStan should have exactly the same result as Stan for the same seed.
>
>
> I think RStan should give exactly the same result as Stan for the same seed.

How could it not?  Is the seed somehow not getting passed?

I think the issue is that rstan resets with a new seed (from R) in lines 499 -- 500

// reset the seed

base_rng.seed(random_seed);

Does it have something to do with the difference
in the way n_eff is calculated?  Coda uses some
spectral thing and RStan uses the split R-hat var+
estimator.

It would be good to have a way to construct a stanfit object from samples*.csv files. 

 

Can you see if the samples themselves are the same?

They are not the same, and also the output from command-line Stan also looks bad. I did this at the command line

goodrich@CYBERPOWERPC:/tmp/SU/Shige$ ./ch7_2 --seed=12345 --chain_id=1 --refresh=-1 --data=Shige.dump

goodrich@CYBERPOWERPC:/tmp/SU/Shige$ mv samples.csv samples_1.csv 

goodrich@CYBERPOWERPC:/tmp/SU/Shige$ ./ch7_2 --chain_id=2 --refresh=-1 --data=Shige.dump

goodrich@CYBERPOWERPC:/tmp/SU/Shige$ mv samples.csv samples_2.csv

goodrich@CYBERPOWERPC:/tmp/SU/Shige$ ./ch7_2 --chain_id=3 --refresh=-1 --data=Shige.dump

goodrich@CYBERPOWERPC:/tmp/SU/Shige$ mv samples.csv samples_3.csv

goodrich@CYBERPOWERPC:/tmp/SU/Shige$ ./ch7_2 --chain_id=4 --refresh=-1 --data=Shige.dump

goodrich@CYBERPOWERPC:/tmp/SU/Shige$ mv samples.csv samples_4.csv

and then this in R

> CLS <- lapply(1:4, FUN = function(chain) read.csv(file=paste0("samples_", chain, ".csv"), comment.char="#"))

> source("Shige.R") # produces the fit object

> options(scipen = 5)

> sapply(CLS, FUN = function(x) colMeans(x[,c(4:8,1)]))

                [,1]         [,2]          [,3]           [,4]

alpha1   105.2152530   105.221594    -0.1216296   105.22197200

alpha2    -0.2451326    71.438486    77.4862292    77.51356160

sigma1     2.7697423     2.776654   105.0612260     2.77333586

sigma2    77.7610272     8.987968     2.7197145     2.68328607

dif      -74.9912841    -6.211314   102.3414980     0.09004965

lp__   -2859.8060600 -1711.659570 -3700.3719800 -1520.92927000

> t(apply(extract(fit, inc_warmup=FALSE), 2:3, mean))

               [,1]           [,2]           [,3]         [,4]

alpha1   105.217908   105.21897243   105.20592224   105.218899

alpha2    59.345505    77.49685872    77.48935035    73.223678

sigma1     2.772667     2.77103205     2.77128175     2.772860

sigma2    19.801809     2.79918901     2.73861566     6.785567

dif      -17.029143    -0.02815696     0.03266609    -4.012706

lp__   -1917.642846 -1524.23329053 -1522.07697366 -1628.017967

Ben

[Bob Carpenter]

On 10/18/12 4:05 PM, Ben Goodrich wrote:
> On Thursday, October 18, 2012 1:18:43 PM UTC-4, Bob Carpenter wrote:
>
> On 10/18/12 12:17 PM, Ben Goodrich wrote:
> > On Thursday, October 18, 2012 11:30:20 AM UTC-4, Jiqiang Guo wrote:

...

> > I think RStan should give exactly the same result as Stan for the same seed.
>
> How could it not? Is the seed somehow not getting passed?
>
>
> I think the issue is that rstan resets with a new seed (from R) in lines 499 -- 500

> |
> // reset the seed
> base_rng.seed(random_seed);

That would certainly do it.

I believe that's what Jiqiang said he'd fixed.

> Does it have something to do with the difference
> in the way n_eff is calculated? Coda uses some
> spectral thing and RStan uses the split R-hat var+
> estimator.
>
> It would be good to have a way to construct a stanfit object from samples*.csv files.

Indeed. I'd like that for testing.

Also, it would be good to have a command-line tool that
runs multiple chains and produces output similar to what
you get with print(fit) in RStan.

> Can you see if the samples themselves are the same?
>
> They are not the same, and also the output from command-line Stan also looks bad. I did this at the command line

OK -- they won't be without the same seed!

> goodrich@CYBERPOWERPC:/tmp/SU/Shige$ ./ch7_2 --seed=12345 --chain_id=1 --refresh=-1 --data=Shige.dump
> goodrich@CYBERPOWERPC:/tmp/SU/Shige$ mv samples.csv samples_1.csv
> goodrich@CYBERPOWERPC:/tmp/SU/Shige$ ./ch7_2 --chain_id=2 --refresh=-1 --data=Shige.dump
> goodrich@CYBERPOWERPC:/tmp/SU/Shige$ mv samples.csv samples_2.csv
> goodrich@CYBERPOWERPC:/tmp/SU/Shige$ ./ch7_2 --chain_id=3 --refresh=-1 --data=Shige.dump
> goodrich@CYBERPOWERPC:/tmp/SU/Shige$ mv samples.csv samples_3.csv
> goodrich@CYBERPOWERPC:/tmp/SU/Shige$ ./ch7_2 --chain_id=4 --refresh=-1 --data=Shige.dump
> goodrich@CYBERPOWERPC:/tmp/SU/Shige$ mv samples.csv samples_4.csv

There's an option for where to dump the samples -- samples.csv
is just the default.

> and then this in R
>
> |
> > CLS <- lapply(1:4, FUN = function(chain) read.csv(file=paste0("samples_", chain, ".csv"), comment.char="#"))
> > source("Shige.R") # produces the fit object
>
> > options(scipen = 5)
> > sapply(CLS, FUN = function(x) colMeans(x[,c(4:8,1)]))
> [,1] [,2] [,3] [,4]
> alpha1 105.2152530 105.221594 -0.1216296 105.22197200
> alpha2 -0.2451326 71.438486 77.4862292 77.51356160
> sigma1 2.7697423 2.776654 105.0612260 2.77333586
> sigma2 77.7610272 8.987968 2.7197145 2.68328607
> dif -74.9912841 -6.211314 102.3414980 0.09004965
> lp__ -2859.8060600 -1711.659570 -3700.3719800 -1520.92927000
> > t(apply(extract(fit, inc_warmup=FALSE), 2:3, mean))
> [,1] [,2] [,3] [,4]
> alpha1 105.217908 105.21897243 105.20592224 105.218899
> alpha2 59.345505 77.49685872 77.48935035 73.223678
> sigma1 2.772667 2.77103205 2.77128175 2.772860
> sigma2 19.801809 2.79918901 2.73861566 6.785567
> dif -17.029143 -0.02815696 0.03266609 -4.012706
> lp__ -1917.642846 -1524.23329053 -1522.07697366 -1628.017967

Assuming the lp__ is correct, the "good" samples are
command-line 4 and RStan 2 and 3? Are those params
sensible?

Is it getting stuck in modes or just not warming up
fast enough?

- Bob

[Ben Goodrich]

On Thursday, October 18, 2012 4:26:41 PM UTC-4, Bob Carpenter wrote:

On 10/18/12 4:05 PM, Ben Goodrich wrote:
> On Thursday, October 18, 2012 1:18:43 PM UTC-4, Bob Carpenter wrote:
>
>     On 10/18/12 12:17 PM, Ben Goodrich wrote:
>      > On Thursday, October 18, 2012 11:30:20 AM UTC-4, Jiqiang Guo wrote:
...
>      > I think RStan should give exactly the same result as Stan for the same seed.
>
>     How could it not?  Is the seed somehow not getting passed?
>
>
> I think the issue is that rstan resets with a new seed (from R) in lines 499 -- 500

> |
> // reset the seed
> base_rng.seed(random_seed);

That would certainly do it.

I believe that's what Jiqiang said he'd fixed.

If so, it hasn't been pushed. 

Assuming the lp__ is correct, the "good" samples are
command-line 4 and RStan 2 and 3?   Are those params
sensible?

The correct parameter values are given in the R file

set.seed(12345)

n1 <- 600                # Number of females

n2 <- 400                # Number of males

alpha1 <- 105                # Population mean of females

alpha2 <- 77.5                # Population mean of males

sigma <- 2.75                # Average population SD of both

n <- n1+n2                # Total sample size

y1 <- rnorm(n1, alpha1, sigma)        # Data for females

y2 <- rnorm(n2, alpha2, sigma)        # Date for males

I originally forgot to specify --warmup=1000 and --iter=21000 in command-line Stan. But the new results are still weird in both cases

>  sapply(CLS, FUN = function(x) colMeans(x[,c(4:8,1)])) 

                [,1]         [,2]          [,3]         [,4]

alpha1   105.2170710     3.878109    -0.6265813   102.266241

alpha2    -0.4422086    77.491623     1.2363986    77.484879

sigma1     2.7687514    99.062619   105.9107940     5.813496

sigma2    78.1094912     2.718003    76.3153247     2.721343

dif      -75.3407397    96.344619    29.5954890     3.092153

lp__   -2860.8249300 -3678.091020 -5034.9440600 -1638.677820

> t(apply(extract(fit, inc_warmup=FALSE), 2:3, mean))

               [,1]           [,2]           [,3]         [,4]

alpha1   105.217908   105.21897243   105.20592224   105.218899

alpha2    59.345505    77.49685872    77.48935035    73.223678

sigma1     2.772667     2.77103205     2.77128175     2.772860

sigma2    19.801809     2.79918901     2.73861566     6.785567

dif      -17.029143    -0.02815696     0.03266609    -4.012706

lp__   -1917.642846 -1524.23329053 -1522.07697366 -1628.017967

Is it getting stuck in modes or just not warming up
fast enough?

I would guess none of the above.

Ben

[Jiqiang Guo]

As I said, if the line of resetting the seed is commented out, RStan and Stan would be the same. 

I have that line of for the purpose it would be reproducible if we take the inits from a previous run and feed the same seed. 

So I would propose to add that line in Stan so that we would have the same results in Stan and RStan.  Of course, I can be convinced to get rid of that line. 

And comparing the results from Stan and RStan if that line is commented out, we need to be careful that the data dumped to Stan 

are different from that directly passed to Stan in RStan as R dump would round the data.  

Ben

--
 
 

[Ben Goodrich]

On Thursday, October 18, 2012 6:01:43 PM UTC-4, Jiqiang Guo wrote:

As I said, if the line of resetting the seed is commented out, RStan and Stan would be the same. 

I have that line of for the purpose it would be reproducible if we take the inits from a previous run and feed the same seed. 

So I would propose to add that line in Stan so that we would have the same results in Stan and RStan.  Of course, I can be convinced to get rid of that line. 

I agree they should be harmonized, but I would be open to which one should change. But the results are different even for the first chain (I used stan_rdump).

[Jiqiang Guo]

Do you use the dumped data in RStan as well?

Jiqiang 

--
 
 

[Bob Carpenter]

Functionally, every time a user supplies the same seed,
they should get the same output (in the same version of
Stan, at least).

Why would we ever RE-set the seed to random if it's
already been set?

- Bob

On 10/18/12 6:01 PM, Jiqiang Guo wrote:
> As I said, if the line of resetting the seed is commented out, RStan and Stan would be the same.
> I have that line of for the purpose it would be reproducible if we take the inits from a previous run and feed the same
> seed.
> So I would propose to add that line in Stan so that we would have the same results in Stan and RStan. Of course, I can
> be convinced to get rid of that line.
>
> And comparing the results from Stan and RStan if that line is commented out, we need to be careful that the data dumped
> to Stan
> are different from that directly passed to Stan in RStan as R dump would round the data.
>

> --
>
>

[Bob Carpenter]

Does the R code calculate the seed for subsequent chains
the same way as the command line? It should be, as that's
the L'Ecuyer "parallel" seeding.

But as Jiqiang points out, the precision can be different
based on the number of decimal places in the dump.

- Bob

> --
>
>

[Ben Goodrich]

Sourcing the attached file is as comparable as I can think of (except for the results).

Ben

[Ben Goodrich]

On Thursday, October 18, 2012 6:41:06 PM UTC-4, Bob Carpenter wrote:

Does the R code calculate the seed for subsequent chains
the same way as the command line?  It should be, as that's
the L'Ecuyer "parallel" seeding.

No, the default PRNG is Mersenne Twister. The L'Ecuyer (1999) algorithm is only an option. And even if they had that option chosen, the seed in stan() defaults to sample.int(.Machine$integer.max) rather than taking the seed from the state of the L'Ecuyer (1999) algorithm. (Which reminds me that I am still skeptical that the L'Ecuyer (1988) algorithm we are using from Boost has the properties of the L'Ecuyer (1999) algorithm).

 

Ben

[Jiqiang Guo]

But you passed the original data list to stan(), not the data sourced.

I am assuming you commented that line in stan_fit.hpp.

I convinced myself that line can be commented out. I will push it soon.

Jiqiang 

--
 
 

<Shige2.txt>

[Ben Goodrich]

Oops. I meant to write

stan_rdump(c("n1", "n2", "y1", "y2"), file = "Shige.dump")

source("Shige.dump")

ch7_2_data <- list(n1=n1, n2=n2, x=x, y1=y1, y2=y2)

I hadn't commented out line 500 of stan_fit.hpp, but when I fix the script as above, the results aren't even the same on the first chain yet.

Ben

[Jiqiang Guo]

The seed is just an integer,I mean I donor care much how it is generated by which algorithm.  In order to let all chains have the same seed, I generated one in R. All the other RNG for sampling  is the same as Stan since I coped over the code.

Jiqiang

--
 
 

[Jiqiang Guo]

Ok. I will demonstrate they are the same.  Give me 20 mins or maybe a bit longer.

Jiqiang 

--
 
 

[Ben Goodrich]

I propose that the default seed for Stan be -1 (which would indicate to generate a seed from the clock time and chain_id) and the default seed for stan() should also be -1. However, I was surprised to find that specifying --seed=-1 works currently

goodrich@CYBERPOWERPC:/tmp/SU/Shige$ ./ch7_2 --data=Shige.dump --seed=-1

STAN SAMPLING COMMAND

data = Shige.dump

init = random initialization

init tries = 1

samples = samples.csv

append_samples = 0

save_warmup = 0

seed = 4294967295 (user specified)

Ben

[Ben Goodrich]

Also, specifying refresh=-1 in command-line Stan seems not to do anything; it looks like it is read as an unsigned int.

goodrich@CYBERPOWERPC:/tmp/SU/Shige$ ./ch7_2 --data=Shige.dump --seed=-1 --refresh=-1

STAN SAMPLING COMMAND

data = Shige.dump

init = random initialization

init tries = 1

samples = samples.csv

append_samples = 0

save_warmup = 0

seed = 4294967295 (user specified)

chain_id = 1 (default)

iter = 2000

warmup = 1000

thin = 1 (default)

equal_step_sizes = 0

leapfrog_steps = -1

max_treedepth = 10

epsilon = -1

epsilon_pm = 0

delta = 0.5

gamma = 0.05

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

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

[Bob Carpenter]

Stan itself uses L'Ecuyer.

Yes, it has the same properties that it's
easy to calculate skip-ahead and the overall
cycle's big enough. I could've gotten it wrong,
of course, but it's pretty straightforward.

- Bob

> --
>
>

[Bob Carpenter]

I'm getting more confused by the message.

You need to replicate the behavior in the
stanc-generated command that moves ahead the
RNG for the chains with IDs other than 0.

- Bob

On 10/18/12 7:37 PM, Jiqiang Guo wrote:
> The seed is just an integer,I mean I donor care much how it is generated by which algorithm. In order to let all chains
> have the same seed, I generated one in R. All the other RNG for sampling is the same as Stan since I coped over the code.
>
> Jiqiang
>

> On Oct 18, 2012, at 7:28 PM, Ben Goodrich <goodri...@gmail.com <mailto:goodri...@gmail.com>> wrote:
>
>> On Thursday, October 18, 2012 6:41:06 PM UTC-4, Bob Carpenter wrote:
>>
>> Does the R code calculate the seed for subsequent chains
>> the same way as the command line? It should be, as that's
>> the L'Ecuyer "parallel" seeding.
>>
>>
>> No, the default PRNG is Mersenne Twister. The L'Ecuyer (1999) algorithm is only an option. And even if they had that
>> option chosen, the seed in stan() defaults to sample.int(.Machine$integer.max) rather than taking the seed from the
>> state of the L'Ecuyer (1999) algorithm. (Which reminds me that I am still skeptical that the L'Ecuyer (1988) algorithm
>> we are using from Boost has the properties of the L'Ecuyer (1999) algorithm).
>> Ben
>>
>> --
>>
>>

> --
>
>

[Bob Carpenter]

On 10/18/12 7:48 PM, Ben Goodrich wrote:
> I propose that the default seed for Stan be -1 (which would indicate to generate a seed from the clock time and
> chain_id) and the default seed for stan() should also be -1. However, I was surprised to find that specifying --seed=-1
> works currently

I think the default behavior should be to choose
a random seed, not use -1.

Seeds are unsigned int in Stan (because they're
uint32_t in Boost).

See:

http://www.boost.org/doc/libs/1_51_0/doc/html/boost/random/ecuyer1988.html

I don't think it matters how the random seed's
generated -- from R, from C++ according to the time,
whatever -- as long as it's available later to recreate.

From stan/gm/command.hpp:

unsigned int random_seed = 0;

...

rng_t base_rng(random_seed);
// (2**50 = 1T samples, 1000 chains)
static boost::uintmax_t DISCARD_STRIDE = static_cast<boost::uintmax_t>(1) << 50;
// DISCARD_STRIDE <<= 50;
base_rng.discard(DISCARD_STRIDE * (chain_id - 1));


typedef boost::ecuyer1988 rng_t;
rng_t base_rng(random_seed);
// (2**50 = 1T samples, 1000 chains)
static boost::uintmax_t DISCARD_STRIDE = static_cast<boost::uintmax_t>(1) << 50;
// DISCARD_STRIDE <<= 50;
base_rng.discard(DISCARD_STRIDE * (chain_id - 1));

The discard skips ahead in the RNG sequence by the
specified number of draws.

- Bob

>
> |
> goodrich@CYBERPOWERPC:/tmp/SU/Shige$ ./ch7_2 --data=Shige.dump --seed=-1
> STAN SAMPLING COMMAND
> data = Shige.dump
> init = random initialization
> init tries = 1
> samples = samples.csv
> append_samples = 0
> save_warmup = 0
> seed = 4294967295 (user specified)
>
> |
>
> Ben
>
> On Thursday, October 18, 2012 7:37:49 PM UTC-4, Jiqiang Guo wrote:
>
> The seed is just an integer,I mean I donor care much how it is generated by which algorithm. In order to let all
> chains have the same seed, I generated one in R. All the other RNG for sampling is the same as Stan since I coped
> over the code.
>
> Jiqiang
>

> On Oct 18, 2012, at 7:28 PM, Ben Goodrich <goodri...@gmail.com <javascript:>> wrote:
>
>> On Thursday, October 18, 2012 6:41:06 PM UTC-4, Bob Carpenter wrote:
>>
>> Does the R code calculate the seed for subsequent chains
>> the same way as the command line? It should be, as that's
>> the L'Ecuyer "parallel" seeding.
>>
>>
>> No, the default PRNG is Mersenne Twister. The L'Ecuyer (1999) algorithm is only an option. And even if they had

>> that option chosen, the seed in stan() defaults to sample.int <http://sample.int>(.Machine$integer.max) rather

>> than taking the seed from the state of the L'Ecuyer (1999) algorithm. (Which reminds me that I am still skeptical
>> that the L'Ecuyer (1988) algorithm we are using from Boost has the properties of the L'Ecuyer (1999) algorithm).
>> Ben
>>
>> --
>>
>>

> --
>
>

[Bob Carpenter]

It's just treated as an int. Nothing in the doc
says it should do anyhting for -1. We should
be throwing an error, I think.

What do you think it should do?

I like the idea of having 0 turn it off,
but haven't implemented it yet.

- Bob

>> had that option chosen, the seed in stan() defaults to sample.int <http://sample.int>(.Machine$integer.max)

>> rather than taking the seed from the state of the L'Ecuyer (1999) algorithm. (Which reminds me that I am still
>> skeptical that the L'Ecuyer (1988) algorithm we are using from Boost has the properties of the L'Ecuyer (1999)
>> algorithm).
>> Ben
>>
>> --
>>
>>

> --
>
>

[Ben Goodrich]

On Thursday, October 18, 2012 8:06:51 PM UTC-4, Bob Carpenter wrote:

It's just treated as an int.  Nothing in the doc
says it should do anyhting for -1.  We should
be throwing an error, I think.

What do you think it should do?

I like the idea of having 0 turn it off,
but haven't implemented it yet.

Well, rstan turns it off for refresh <= 0, so it would probably be good to be consistent. And I can't think of another reason why someone would supply a negative number.

Ben

 

[Jiqiang Guo]

Well, for the time being, I failed.  But I think in the morning, I made it happened.  I am working on it and will make sure they would be the same though. 

In addition, when iter=21000, the thin in Stan is not 1. 

--

Jiqiang

[Jiqiang Guo]

Trying to make what I said more clear. 

I am just saying all the chains in RStan are using the same seed but different chain_id.  That is, in the following code, random_seed is the same for all chains; but chain_id's are different, the default is from 1 to $chains. 

      typedef boost::ecuyer1988 rng_t;
      rng_t base_rng(random_seed);
      // (2**50 = 1T samples, 1000 chains)
      static boost::uintmax_t DISCARD_STRIDE = static_cast<boost::uintmax_t>(1) << 50;
      // DISCARD_STRIDE <<= 50;
      base_rng.discard(DISCARD_STRIDE * (chain_id - 1));

--

[Ben Goodrich]

On Thursday, October 18, 2012 8:04:43 PM UTC-4, Bob Carpenter wrote:

On 10/18/12 7:48 PM, Ben Goodrich wrote:
> I propose that the default seed for Stan be -1 (which would indicate to generate a seed from the clock time and
> chain_id) and the default seed for stan() should also be -1. However, I was surprised to find that specifying --seed=-1
> works currently

I think the default behavior should be to choose
a random seed, not use -1.

Seeds are unsigned int in Stan (because they're
uint32_t in Boost).

If stan() is going to be passing a seed, then stan() is limited to signed integers because R doesn't have an unsigned integer type. So, there might be sequences that someone could do in Stan that wouldn't be possible to replicate with rstan. Maybe reserve 0 as the code for "generate a seed from the time and chain_id"?

Ben

[Ben Goodrich]

I've not been able to find anything that comes out and enthusiastically says algorithms like his 1988 algorithm are okay for statistical computing. L'Ecuyer's most recent survey

http://scholar.google.com/scholar?cluster=16324058590511781724&hl=en&as_sdt=0,33

says that the period should be much bigger than the square of the number of pseudo-random numbers you need. I guess that depends on the number of parameters in our case. Is it conceivable that we would need 2 billion? In the last section, he says "Computer implementations that this author can suggest for the moment include ... " a few of his from the 1990s, the Mersenne-Twister ones, and a Fortran one.

Still don't really know what to make of all this.

Ben

[Jiqiang Guo]

I do not know what is going on now.  I guess this morning, I just looked at samples before warmup. 

There might be a big issue in RStan.  It seems that the NUTS adaptation in RStan has never been turned off though I did call sampler.adapt_off() to turn if off after warmup.  I looked at gm/command.hpp and rstan/stanfit.hpp many times, but could not figure out where I did it wrong. 

--

Jiqiang

[Jiqiang Guo]

Ok. Finally I found the bug in RStan.  The bug/typo is that the value of epsilon_pm was -1 (the default value of epsilon) instead of 0.  I do not know how bad that would influence the sampling.  I guess we need to release a new version of RStan soon.

I'll push the code soon. 

--

Jiqiang

[Jiqiang Guo]

I changed the original 'Shige2.txt' a bit so all the results should be exactly the same in RStan and Stan. 

The final results for running the script is because samples printed out by Stan have less digits after decimal points. 

--

Jiqiang

[Ben Goodrich]

On Friday, October 19, 2012 10:38:56 PM UTC-4, Jiqiang Guo wrote:

I changed the original 'Shige2.txt' a bit so all the results should be exactly the same in RStan and Stan. 

The final results for running the script is because samples printed out by Stan have less digits after decimal points. 

Great! Now we are down to the problem that sigma2 in the first chain is equally terrible in both cases. That and the salm model either goes badly or with seed = 12345 yields

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

But that is equally true for rstan and Stan. Replication script attached. You will need to modify STAN_HOME for it to (not) work.

Ben

[Bob Carpenter]

Let's target a release next week.

I'm most of the way through adding upper-/lower-bound constraints
for vectors and matrices, but even if I don't
finish that, we have a bunch of issues that make
it worth releasing 1.0.3 sooner rather than later.

- Bob

On 10/19/12 10:20 PM, Jiqiang Guo wrote:
> Ok. Finally I found the bug in RStan. The bug/typo is that the value of epsilon_pm was -1 (the default value of
> epsilon) instead of 0. I do not know how bad that would influence the sampling. I guess we need to release a new
> version of RStan soon.
>
> I'll push the code soon.
>
> --
> Jiqiang
>

> On Thu, Oct 18, 2012 at 11:27 PM, Jiqiang Guo <guo...@gmail.com <mailto:guo...@gmail.com>> wrote:
>
> I do not know what is going on now. I guess this morning, I just looked at samples before warmup.
>
> There might be a big issue in RStan. It seems that the NUTS adaptation in RStan has never been turned off though I
> did call sampler.adapt_off() to turn if off after warmup. I looked at gm/command.hpp and rstan/stanfit.hpp many
> times, but could not figure out where I did it wrong.
>
> --
> Jiqiang
>
>
> On Thu, Oct 18, 2012 at 8:40 PM, Jiqiang Guo <guo...@gmail.com <mailto:guo...@gmail.com>> wrote:
>
> Well, for the time being, I failed. But I think in the morning, I made it happened. I am working on it and
> will make sure they would be the same though.
>
> In addition, when iter=21000, the thin in Stan is not 1.
>
> --
> Jiqiang
>

> On Thu, Oct 18, 2012 at 7:41 PM, Jiqiang Guo <guo...@gmail.com <mailto:guo...@gmail.com>> wrote:
>
> Ok. I will demonstrate they are the same. Give me 20 mins or maybe a bit longer.
>
> Jiqiang
>

> --
>
>

[Bob Carpenter]

It's very configurable as long as we can implement
the Boost concept of a base random number generator.

The question I see is whether we need an RNG for
which it's efficient to advance the position in the
sequence. This is what the one we're using lets you
do and I thought you (Ben) thought this was important.

So the issue would be how to generate a seed for
each chain using the Mersenne Twister or other RNG.

Is there some way to answer this question empirically?

Without forward simulation in generated quantities (not
implemented yet), we're bounded by:

(dimensions + 1) * iterations * chains

Still, 2G isn't that much headroom. 50K iterations of
20K parameters would do it in a single chain.

But ecuyer1988 we're using is 2^61.

The citation is:

Efficient and Portable Combined Random Number Generators", Pierre L'Ecuyer, Communications of the ACM, Vol. 31, No. 6,
June 1988, pp. 742-749, 774

(I won't be able to read the paywalled ref until I
get back to Columbia.)

- Bob

On 10/18/12 10:14 PM, Ben Goodrich wrote:
> I've not been able to find anything that comes out and enthusiastically says algorithms like his 1988 algorithm are okay
> for statistical computing. L'Ecuyer's most recent survey
>
> http://scholar.google.com/scholar?cluster=16324058590511781724&hl=en&as_sdt=0,33
>
> says that the period should be much bigger than the square of the number of pseudo-random numbers you need. I guess that
> depends on the number of parameters in our case. Is it conceivable that we would need 2 billion? In the last section, he
> says "Computer implementations that this author can suggest for the moment include ... " a few of his from the 1990s,
> the Mersenne-Twister ones, and a Fortran one.
>
> Still don't really know what to make of all this.
>
> Ben
>
> On Thursday, October 18, 2012 7:56:30 PM UTC-4, Bob Carpenter wrote:
>
> Stan itself uses L'Ecuyer.
>
> Yes, it has the same properties that it's
> easy to calculate skip-ahead and the overall
> cycle's big enough. I could've gotten it wrong,
> of course, but it's pretty straightforward.
>
> - Bob
>
>
> On 10/18/12 7:28 PM, Ben Goodrich wrote:
> > On Thursday, October 18, 2012 6:41:06 PM UTC-4, Bob Carpenter wrote:
> >
> > Does the R code calculate the seed for subsequent chains
> > the same way as the command line? It should be, as that's
> > the L'Ecuyer "parallel" seeding.
> >
> >
> > No, the default PRNG is Mersenne Twister. The L'Ecuyer (1999) algorithm is only an option. And even if they had that

> > option chosen, the seed in stan() defaults to sample.int <http://sample.int>(.Machine$integer.max) rather than

> taking the seed from the
> > state of the L'Ecuyer (1999) algorithm. (Which reminds me that I am still skeptical that the L'Ecuyer (1988)
> algorithm
> > we are using from Boost has the properties of the L'Ecuyer (1999) algorithm).
> > Ben
> >
> > --
> >
> >
>

> --
>
>

[Bob Carpenter]

They're both 32 bits on all the platforms we currently
care about and there's a 1-to-1 relationship defined
by just reinterpreting the bits :-)

To see what a signed int looks like as unsigned:

int x;
std::cout << x << "=" << static_cast<unsigned int>(x);

and similarly the other way around.

- Bob

[Ben Goodrich]

On Saturday, October 20, 2012 4:44:55 PM UTC-4, Bob Carpenter wrote:

Without forward simulation in generated quantities (not
implemented yet), we're bounded by:

   (dimensions + 1) * iterations * chains

Does each leapfrog step not require any randomness?

 

(I won't be able to read the paywalled ref until I
get back to Columbia.)

Not that you need to read anything now, now, now. But if you goto

http://scholar.google.com/

sign in with your gmail account, click settings in the top right, and then library links in the top left, it will allow you to enter Columbia University. After you click save, results on Google Scholar will have an "e-Link @ Columbia" which will allow you to access electronic articles with your Columbia ID and password.

Ben

[Ben Goodrich]

On Friday, October 19, 2012 11:39:23 PM UTC-4, Ben Goodrich wrote:

On Friday, October 19, 2012 10:38:56 PM UTC-4, Jiqiang Guo wrote:

I changed the original 'Shige2.txt' a bit so all the results should be exactly the same in RStan and Stan. 

The final results for running the script is because samples printed out by Stan have less digits after decimal points. 

Great! Now we are down to the problem that sigma2 in the first chain is equally terrible in both cases.

In this case, the means are 105 and 77.5, so starting either at 0 or 0 in expectation is a long way from the truth. It works if you standardize the input. So, maybe this is not a blocker (haven't looked into the salm problem yet though), but it would be nice if Stan could adapt a little better.

Ben

 

[Marcus Brubaker]

That's brilliant!  Thanks for sharing :)

Cheers,
Marcus

[Ben Goodrich]

On Saturday, October 20, 2012 4:44:55 PM UTC-4, Bob Carpenter wrote:

It's very configurable as long as we can implement
the Boost concept of a base random number generator.

The question I see is whether we need an RNG for
which it's efficient to advance the position in the
sequence.  This is what the one we're using lets you
do and I thought you (Ben) thought this was important.

That is important. I mentioned in another thread that

Because I don't understand this stuff well, I would feel better if we used something from 
http://trng.berlios.de/
as the PRNG. It is BSD, it is small, it is designed to be like Boost, it formed the basis for a proposal to get it into the C++11 standard, and the designers have explicitly thought about parallelism. So, to me, it is worth another dependency.

 

Their PRNGs all have methods to split the stream and jump ahead:

void split(unsigned int, unsigned int);

void jump2(unsigned int);

void jump(unsigned long long);

It also has inverse CDFs.

Ben

[Bob Carpenter]

On 10/20/12 6:24 PM, Marcus Brubaker wrote:

> On Sat, Oct 20, 2012 at 5:53 PM, Ben Goodrich <goodri...@gmail.com <mailto:goodri...@gmail.com>> wrote:
>
> On Saturday, October 20, 2012 4:44:55 PM UTC-4, Bob Carpenter wrote:
>
> Without forward simulation in generated quantities (not
> implemented yet), we're bounded by:
>
> (dimensions + 1) * iterations * chains
>
>
> Does each leapfrog step not require any randomness?
>
> (I won't be able to read the paywalled ref until I
> get back to Columbia.)

No, there's just the initial randomization of the momentum
vector. Then it's deterministic until the Metropolis step,
which may or may not require RNG in [0,1].

When we do forward simulation in generated quantities,
there may be additional randomness per iteration (not per
leapfrog step).

>
> Not that you need to read anything now, now, now. But if you goto
>
> http://scholar.google.com/
>
> sign in with your gmail account, click settings in the top right, and then library links in the top left, it will
> allow you to enter Columbia University. After you click save, results on Google Scholar will have an "e-Link @
> Columbia" which will allow you to access electronic articles with your Columbia ID and password.
>
>
> That's brilliant! Thanks for sharing :)

Indeed. Thanks.

- Bob

[Bob Carpenter]

On 10/20/12 7:01 PM, Ben Goodrich wrote:
> On Saturday, October 20, 2012 4:44:55 PM UTC-4, Bob Carpenter wrote:
>
> It's very configurable as long as we can implement
> the Boost concept of a base random number generator.
>
> The question I see is whether we need an RNG for
> which it's efficient to advance the position in the
> sequence. This is what the one we're using lets you
> do and I thought you (Ben) thought this was important.
>
>
> That is important. I mentioned in another thread

I recall that.

> that
>
> Because I don't understand this stuff well,

Speaking of which, I did go with one that jumped ahead
to support parallelism. I just didn't realize there
were multiple L'Ecuyer randomizers with that property!

> ... I would feel better if we used something from

> http://trng.berlios.de/
> as the PRNG. It is BSD, it is small, it is designed to be like Boost, it formed the basis for a proposal to get it
> into the C++11 standard, and the designers have explicitly thought about parallelism. So, to me, it is worth another
> dependency.
>
> Their PRNGs all have methods to split the stream and jump ahead:
>
> void split(unsigned int, unsigned int);
> void jump2(unsigned int);
> void jump(unsigned long long);

It's OK with me.

The biggest downside is that it adds another dependency
and library to keep up to date.

If it doesn't implement the Boost concept, we could
either (a) adapt it in another class so that it does,
or (b) reimplement the random draws in the mcmc
samplers for normal(0,1) and uniform(0,1), which looks
easy given:

> It also has inverse CDFs.

Assuming that's the best way to implement normal(0,1).

Boost has inverse CDFs, too. I do want to start implementing
forward generation, and we might as well decide on whether
to rip out the current randomizers and replace them before
starting the forward generation.

The main issue for CDFs to support truncation/censoring
is computing gradients.

- Bob

[Jiqiang Guo]

In next release of rstan, the seed could be a string of digits, say, seed = "12345".  As a result, rstan 

is able to use any seed that is used in Stan. 

--

Jiqiang

[Ben Goodrich]

Would it make sense to have init default to 2 and draw initial values from the (-init, init) interval? That would make it easy to get initial values from a wider or narrow interval and still preserve the existing behavior when init is zero.

Ben

[Bob Carpenter]

n 10/23/12 1:32 AM, Ben Goodrich wrote:
> Would it make sense to have init default to 2 and draw initial values from the (-init, init) interval? That would make
> it easy to get initial values from a wider or narrow interval and still preserve the existing behavior when init is zero.

Yes, this would be very easy.

It would also be easy to sample from a distribution
as Marcus suggested at some point.

It's all in the stan/gm/command.hpp code, starting
around line 351:

init_val = "random initialization"; // for I/O
// init_rng generates uniformly from -2 to 2
boost::random::uniform_real_distribution<double>
init_range_distribution(-2.0,2.0);
boost::variate_generator<rng_t&,
boost::random::uniform_real_distribution<double> >
init_rng(base_rng,init_range_distribution);

params_i = std::vector<int>(model.num_params_i(),0);
params_r = std::vector<double>(model.num_params_r());

// retry inits until get a finite log prob value
std::vector<double> init_grad;
static int MAX_INIT_TRIES = 100;
for (num_init_tries = 1; num_init_tries <= MAX_INIT_TRIES; ++num_init_tries) {
for (size_t i = 0; i < params_r.size(); ++i)
params_r[i] = init_rng();

All you need to do is extract the init value from the command
line and give it a default, then use it instead of 2.0 in
the above. Using a different distribution other than uniform's
also easy -- that's just the init_range_distribution spec.

- Bob