How to speed up my stan code?

[Li Pei Wang]

hello,

The following paragraph is my stan code.

I ran 2-chain with 2-cores and iteration 150000 times.

It took around 69 hrs.

I really want to know how can I speeding up my code.

Thank you.

 data

{

int <lower=0> n;//206

int <lower=0> N;//2695

cov_matrix[3] S;

vector<lower=0>[n] stage;//early stage and advanced stage

int<lower=0> ID[N];

vector[N] observetime;

vector[N] log10ca125;

}

parameters

{

real u1; // sigma

real<lower=0> u2; // sigma_stage

real u11; // beta1

real<lower=0> u12; // beta1_stage

real beta2; //beta2

real beta3; //beta3

real delta; //skewness parameters of random error

vector[3] deltab; //skewness parameters of random effect

matrix<lower=0>[n,3]  wb; //wb = |X_0b|~|N(0,1)|

cov_matrix[3] cov;

matrix[n,3] b; //random effect 

real<lower=3> nu; //degree of freedom

vector<lower=0>[n] ue; //latent random variables generated by ST distribution 

vector<lower=0>[N] w;

vector[N] Y_pred;

}

transformed parameters

{

vector[n] beta1; // Stage info for extra random effects mean

vector<lower=0>[n] sigma; // Standard deviation

matrix[n,3] mub;//mean of random effect matrix

real d1;

real<lower=0> mue;

//----------------------------------------

vector[n] p1;

vector[n] p2;

vector[n] p3;

//----------------------------------------

vector[N] dm;

vector[N] mu;

vector[N] aau;

vector[N] res;

vector[N] sue;

real sten;

//----------------------------------------

beta1 = u11 + u12 * stage;

sigma = exp(u1 + u2 * stage);

for(i in 1:3)

{

mub[,i] = deltab[i] * (wb[,i] - sqrt(2/3.14159));

}

//means of error term (ST distribution)

mue = exp(lgamma(0.5*(nu-1)) -lgamma(0.5*nu))*sqrt(nu/3.14159);

d1 = nu/2;

sten = sqrt(10);

// Mixture of nonlinear mixed-effects models with two components using ST distribution 

p1 = beta1 + b[,1];

p2 = beta2 + b[,2];

p3 = beta3 + b[,3];

 

sue = sqrt(ue[ID]);

aau = sue .* sigma[ID];

dm = ( p1[ID] + p2[ID] .* observetime + p3[ID] .* observetime .* observetime)/log(10);

mu = dm + delta*(w-mue);

res = log10ca125 - mu;

}

//----------------------------------------

model

{

//degree of freedom

nu ~ exponential(0.1);

// Linear regression

u1 ~ normal(0,sten);

u2 ~ normal(0,sten);

u11 ~ normal(0,sten);

u12 ~ normal(0,sten);

//skewness parameters 

  delta ~ normal(0,sten);  //random error 

  deltab ~ normal(0,sten); //random effect 

//Random effects of response model

for(i in 1:3)

wb[,i] ~ normal(0,1);

//Variance-covariance matrix

  cov ~ inv_wishart(3,S);

//random effect

for(s in 1:n)

b[s,] ~ multi_normal(mub[s,],cov);

//Latent random variables generated by ST distribution 

  ue ~ inv_gamma(d1,d1);

w ~ normal(0,sue);

//Coefficients 

beta2 ~ normal(0,sten);

beta3 ~ normal(0,sten);

//Generate ST random variable y from the SN tandom variable z.

log10ca125 ~ normal(mu, aau);

Y_pred ~ normal(mu, aau);

}

[Bob Carpenter]

First, you probably don't need 150K iterations with NUTS.

Vectorize the multi-normal.

> for(s in 1:n)
> b[s,] ~ multi_normal(mub[s,],cov);

as

b ~ multi_normal(mub, cov);

That means you have to declare b and mub as

vector[3] b[n];

and

vector[3] mub[n];

We recommend a scaled correlation matrix parameterized with
a Cholesky factor---it's much more efficient. We don't yet
have very efficient Wisharts and inverse Wisharts, but we
will might eventually if someone writes the Cholesky-factor
parameter versions.

And vectorize

> for(i in 1:3)
> wb[,i] ~ normal(0,1);

as

to_vector(wb) ~ normal(0, 1);

If you don't care about all those transformed parameters, just
make them local variables in the model. Then they won't
get saved. 150K draws might take a lot of memory if n is large
and you save all those transforms.

- 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.
> To post to this group, send email to stan-...@googlegroups.com.
> For more options, visit https://groups.google.com/d/optout.

[Bob Carpenter]

Oh, and you probably want to use a non-centered parameterization
for the hierarchical model. There's an explanation in the
manual.

- Bob

[Li Pei Wang]

Thank you, Bob!

and I still have question about that:

If you don't care about all those transformed parameters, just
make them local variables in the model.  Then they won't
get saved.  150K draws might take a lot of memory if n is large
and you save all those transforms.

Is it means that if I dont need to save the parameter. Then I can just do something like this

 

model

{

vector[3] mub[n];//mean of random effect matrix

real d1;

real<lower=0> mue;

real sten;

vector[N] dm;

vector[N] mu;

vector[N] aau;

vector[N] sue;

for(i in 1:3)

mub = deltab .* (wb[,i] - sqrt(2/3.14159));

//means of error term (ST distribution)

mue = exp(lgamma(0.5*(nu-1)) -lgamma(0.5*nu))*sqrt(nu/3.14159);

d1 = nu/2;

sten = sqrt(10);

sue = sqrt(ue[ID]);

aau = sue .* sigma[ID];

dm = ( p1[ID] + p2[ID] .* observetime + p3[ID] .* observetime .* observetime)/log(10);

  mu = dm + delta*(w-mue);

//degree of freedom

nu ~ exponential(0.1);

// Linear regression

u1 ~ normal(0,sten);

u2 ~ normal(0,sten);

u11 ~ normal(0,sten);

u12 ~ normal(0,sten);

//skewness parameters 

  delta ~ normal(0,sten);  //random error 

  deltab ~ normal(0,sten); //random effect 

//Random effects of response model

to_vector(wb) ~ normal(0,1);

//Variance-covariance matrix

  cov ~ inv_wishart(3,S);

//random effect

b ~ multi_normal(mub,cov);

[Daniel Lee]

Before talking about speed, did your model converge? And what's the scale of your data?

Daniel

[Li Pei Wang]

hi Daniel,

thanks for your reply.

Actually, after iteration 150000 times, it didnt converge in stan, but converge in WinBUGS.

I am confuse about that.

They both run 150000 times, WinBUGS takes 3 hrs (converge), stan takes 15 hrs (do not converge).

The scale of data is around 0 to 5.

 

Daniel Lee於 2017年2月15日星期三 UTC+8下午9時56分39秒寫道：

[Bob Carpenter]

We don't know of any models where JAGS/BUGS converge but Stan
doesn't, so my first guess is that the models are not the same.

Did you run multiple chains from random starting points in WinBUGS?
I'm not sure what it does by default to test convergence, but Stan's
tests are more stringent by design.

- Bob

[Andrew Gelman]

Nothing should be taking 150000 iterations.  I can't be sure but I'm guessing there are big problems with your model if it's needing so many iterations.  The model seems simple enough that I can't see why it wouldn't converge after 2000 iterations!

Also, just a minor thing:  in Stan we have pi().  And sqrt(2/pi) can be computed in the transformed data block.  But I think there must be some simple problems with your model somewhere (in Bugs or Stan).

A

[Matt Espe]

One issue is Y_pred should not be in the model block - this should be in the generated quantities. 

Matt

[Li Pei Wang]

Many thanks for the prompt reply.

I checked my model. Because that I did something wrong on the data set, I get the wrong result. Thanks for Bob's comment, now my model is correct and same to the WinBUGS model.

I still have a question.

Is it possible that the model didn't converge at some random seed, but converge at some random seed?

The following result is random seed = 10, iteration = 30000, warmup = 15000, thin = 50, Chain = 4.

           mean se_mean   sd  2.5%   25%   50%   75% 97.5% n_eff Rhat

nu         3.33    0.01 0.26  3.01  3.13  3.27  3.46  3.99   479 1.01

beta2     -0.06    0.00 0.01 -0.08 -0.07 -0.06 -0.06 -0.05    41 1.09

beta3      0.00    0.00 0.00  0.00  0.00  0.00  0.00  0.00    24 1.19

u1        -3.48    0.08 0.21 -3.90 -3.62 -3.48 -3.32 -3.08     6 1.34

u2         0.37    0.04 0.20  0.04  0.22  0.36  0.51  0.82    28 1.11

u11        3.25    0.01 0.09  3.08  3.19  3.25  3.31  3.42    40 1.11

u12        0.28    0.01 0.09  0.11  0.22  0.28  0.34  0.46    85 1.02

delta      0.52    0.00 0.02  0.48  0.51  0.52  0.53  0.55    20 1.23

deltab[1]  0.59    0.07 0.27 -0.04  0.42  0.65  0.78  0.97    15 1.22

deltab[2] -0.06    0.00 0.02 -0.10 -0.08 -0.07 -0.05  0.00    26 1.11

deltab[3]  0.00    0.00 0.00  0.00  0.00  0.00  0.00  0.00   533 1.00

cov[1,1]   0.24    0.02 0.10  0.07  0.16  0.23  0.32  0.44    16 1.17

cov[1,2]  -0.01    0.00 0.00 -0.02 -0.01 -0.01  0.00  0.00    70 1.05

cov[1,3]   0.00    0.00 0.00  0.00  0.00  0.00  0.00  0.00   319 1.00

cov[2,1]  -0.01    0.00 0.00 -0.02 -0.01 -0.01  0.00  0.00    70 1.05

cov[2,2]   0.00    0.00 0.00  0.00  0.00  0.00  0.00  0.00    48 1.06

cov[2,3]   0.00    0.00 0.00  0.00  0.00  0.00  0.00  0.00   453 1.01

cov[3,1]   0.00    0.00 0.00  0.00  0.00  0.00  0.00  0.00   319 1.00

cov[3,2]   0.00    0.00 0.00  0.00  0.00  0.00  0.00  0.00   453 1.01

cov[3,3]   0.00    0.00 0.00  0.00  0.00  0.00  0.00  0.00   232 1.01

the result seems getting convergence.

So, I restart the iteration, and change the random seed.

seed = 921, iteration = 60000, warmup = 30000, thin = 50, Chain = 4.

and get the following result.

           mean se_mean   sd  2.5%   25%   50%   75% 97.5% n_eff Rhat

nu         3.36    0.01 0.28  3.01  3.15  3.31  3.52  4.02   363 1.01

beta2     -0.06    0.00 0.01 -0.08 -0.06 -0.06 -0.05 -0.04     6 1.21

beta3      0.00    0.00 0.00  0.00  0.00  0.00  0.00  0.00    14 1.28

u1        -3.84    0.38 0.56 -4.92 -4.24 -3.68 -3.41 -3.03     2 3.25

u2         0.81    0.45 0.67  0.05  0.30  0.58  1.35  2.07     2 3.45

u11        3.27    0.01 0.08  3.11  3.21  3.27  3.32  3.43    47 1.08

u12        0.28    0.01 0.08  0.09  0.22  0.29  0.33  0.44   112 1.07

delta      0.53    0.01 0.02  0.50  0.52  0.53  0.55  0.56     4 1.43

deltab[1]  0.59    0.13 0.43 -0.66  0.51  0.73  0.87  1.06    11 1.34

deltab[2] -0.06    0.00 0.03 -0.10 -0.08 -0.06 -0.04  0.01    42 1.12

deltab[3]  0.00    0.00 0.00  0.00  0.00  0.00  0.00  0.00   333 1.00

cov[1,1]   0.20    0.03 0.10  0.06  0.13  0.18  0.28  0.42    10 1.29

cov[1,2]  -0.01    0.00 0.00 -0.02 -0.01 -0.01  0.00  0.00     8 1.20

cov[1,3]   0.00    0.00 0.00  0.00  0.00  0.00  0.00  0.00   358 1.01

cov[2,1]  -0.01    0.00 0.00 -0.02 -0.01 -0.01  0.00  0.00     8 1.20

cov[2,2]   0.00    0.00 0.00  0.00  0.00  0.00  0.00  0.00    77 1.06

cov[2,3]   0.00    0.00 0.00  0.00  0.00  0.00  0.00  0.00   162 1.03

cov[3,1]   0.00    0.00 0.00  0.00  0.00  0.00  0.00  0.00   358 1.01

cov[3,2]   0.00    0.00 0.00  0.00  0.00  0.00  0.00  0.00   162 1.03

cov[3,3]   0.00    0.00 0.00  0.00  0.00  0.00  0.00  0.00   150 1.03

It seems far from the result of seed = 10. 

So, it's possible that the model will converge at some random seed, but not at all random seed?

Lipei

[Michael Betancourt]

The chains converge or they don’t — if some seeds behave differently

than others then _none_ of the chains converged in the first place.

[Bob Carpenter]

You're not converginghere with Rhat >> 1.05.
If you want to repost the model, we can try to help.

- Bob