model fails when rerun

[Peter Li]

When I run the attached stan model and data once using RStan, I get good results. However, if I run it again, I get terrible results.

Here is the code:

roachdata <- read.csv ("roachdata.csv")
attach(roachdata)

if (!exists("roaches.sm")) {
    if (file.exists("roaches.sm.RData")) {
        load("roaches.sm.RData", verbose = TRUE)
    } else {
        rt <- stanc("roaches.stan", model_name = "roaches")
        roaches.sm <- stan_model(stanc_ret = rt)
        save(roaches.sm, file = "roaches.sm.RData")
    }
}

dataList.1 <- list(N=length(y), y=y,roach1=roach1,treatment=treatment,exposure2=exposure2,senior=senior)
roaches.sf1 <- sampling(roaches.sm, dataList.1)
print(roaches.sf1)

Here is the output:

> roachdata <- read.csv ("roachdata.csv")
> attach(roachdata)
The following object(s) are masked from 'roachdata (position 3)':

    exposure2, roach1, senior, treatment, X, y
The following object(s) are masked from 'roachdata (position 4)':

    exposure2, roach1, senior, treatment, X, y
The following object(s) are masked from 'roachdata (position 5)':

    exposure2, roach1, senior, treatment, X, y
> 
> if (!exists("roaches.sm")) {
+     if (file.exists("roaches.sm.RData")) {
+         load("roaches.sm.RData", verbose = TRUE)
+     } else {
+         rt <- stanc("roaches.stan", model_name = "roaches")
+         roaches.sm <- stan_model(stanc_ret = rt)
+         save(roaches.sm, file = "roaches.sm.RData")
+     }
+ }

TRANSLATING MODEL 'roaches' FROM Stan CODE TO C++ CODE NOW.
COMPILING THE C++ CODE FOR MODEL 'roaches' NOW.

dataList.1 <- list(N=length(y), y=y,roach1=roach1,treatment=treatment,exposure2=exposure2,senior=senior)
roaches.sf1 <- sampling(roaches.sm, dataList.1)
print(roaches.sf1)

> 
> dataList.1 <- list(N=length(y), y=y,roach1=roach1,treatment=treatment,exposure2=exposure2,senior=senior)
> roaches.sf1 <- sampling(roaches.sm, dataList.1)
SAMPLING FOR MODEL 'roaches' NOW (CHAIN 1).
Iteration: 2000 / 2000 [100%]  (Sampling)

SAMPLING FOR MODEL 'roaches' NOW (CHAIN 2).
Iteration: 2000 / 2000 [100%]  (Sampling)

SAMPLING FOR MODEL 'roaches' NOW (CHAIN 3).
Iteration: 2000 / 2000 [100%]  (Sampling)

SAMPLING FOR MODEL 'roaches' NOW (CHAIN 4).
Iteration: 2000 / 2000 [100%]  (Sampling)

> print(roaches.sf1)
Inference for Stan model: roaches.
4 chains, each with iter=2000; warmup=1000; thin=1; 
post-warmup draws per chain=1000, total post-warmup draws=4000.

           mean se_mean  sd    2.5%     25%     50%     75%   97.5% n_eff Rhat
beta[1]     3.1       0 0.0     3.0     3.1     3.1     3.1     3.1  1623    1
beta[2]     0.0       0 0.0     0.0     0.0     0.0     0.0     0.0  1666    1
beta[3]    -0.5       0 0.0    -0.6    -0.5    -0.5    -0.5    -0.5  1901    1
beta[4]    -0.4       0 0.0    -0.4    -0.4    -0.4    -0.4    -0.3  1591    1
lp__    17603.9       0 1.5 17600.1 17603.2 17604.3 17605.0 17605.7  1270    1

Samples were drawn using NUTS2 at Tue Jul 30 14:27:23 2013.
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).
> 
> 
> dataList.1 <- list(N=length(y), y=y,roach1=roach1,treatment=treatment,exposure2=exposure2,senior=senior)
> roaches.sf1 <- sampling(roaches.sm, dataList.1)
SAMPLING FOR MODEL 'roaches' NOW (CHAIN 1).
Iteration:  200 / 2000 [ 10%]  (Adapting)print(roaches.sf1)
Iteration: 2000 / 2000 [100%]  (Sampling)

SAMPLING FOR MODEL 'roaches' NOW (CHAIN 2).
Iteration: 2000 / 2000 [100%]  (Sampling)

SAMPLING FOR MODEL 'roaches' NOW (CHAIN 3).
Iteration: 2000 / 2000 [100%]  (Sampling)

SAMPLING FOR MODEL 'roaches' NOW (CHAIN 4).
Iteration: 2000 / 2000 [100%]  (Sampling)

> print(roaches.sf1)
Inference for Stan model: roaches.
4 chains, each with iter=2000; warmup=1000; thin=1; 
post-warmup draws per chain=1000, total post-warmup draws=4000.

                  mean se_mean  sd          2.5%            25%     50%     75%
beta[1]   2.500000e+00     0.8 1.1   6.00000e-01   2.400000e+00     3.1     3.1
beta[2]   3.000000e-01     0.4 0.6   0.00000e+00   0.000000e+00     0.0     0.3
beta[3]  -9.000000e-01     0.4 0.6  -2.00000e+00  -1.000000e+00    -0.5    -0.5
beta[4]  -7.000000e-01     0.4 0.5  -1.60000e+00  -8.000000e-01    -0.4    -0.4
lp__    -1.061407e+250     Inf Inf -4.24563e+250 -1.061407e+250 17603.8 17604.9
          97.5% n_eff   Rhat
beta[1]     3.1     2   63.0
beta[2]     1.3     2 7813.7
beta[3]    -0.5     2   33.2
beta[4]    -0.3     2   19.0
lp__    17605.7  4000    NaN

Samples were drawn using NUTS2 at Tue Jul 30 14:27:45 2013.
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).

Peter

[Daniel Lee]

Have you looked at the initial values for the second run? If you have the seed laying around, you might want to check the gradients for the model with how it's being initialized.


Daniel

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

[Peter Li]

Ah, I see. It has nothing to do with running the model twice. It has to
do with the initial values. Thanks.


Peter