Stan and cmdstan running slow

Modeling
1

Hi all,

I am using the code below. This code was modified to use ver weak priors since I was having some issues of convergence before. However, this takes too long, since it has been running all night and it did reach 30% in Chain 1. Moreover, by using cmdstanr besides stan_model. Any ideas how can I improve on this? Is this normal that it takes so long?

Chain 1: Gradient evaluation took 0.022591 seconds
Chain 1: 1000 transitions using 10 leapfrog steps per transition would take 225.91 seconds.
Chain 1: Adjust your expectations accordingly!
Chain 1:
Chain 1:
Chain 1: Iteration: 1 / 2000 [ 0%] (Warmup)
Chain 1: Iteration: 200 / 2000 [ 10%] (Warmup)
Chain 1: Iteration: 400 / 2000 [ 20%] (Warmup)
Chain 1: Iteration: 600 / 2000 [ 30%] (Warmup)

functions {
  vector model3(real t, vector y, real pAB, real uAB, real uASC) {
    vector[2] dydt;  
    dydt[1] = pAB*y[2]-uAB*y[1];
    dydt[2] = -uASC*y[2];
    
   return dydt;
  }
}
data {
  int <lower=1> nobs;
  real t0;
  vector[2] y0;
  real ts[nobs];
  int <lower=1> indivs;
  real <lower=0> antib[nobs];
  //real <lower=1, upper=indivs> subj[nobs];
}
parameters {
  real <lower=0> pAB0;
  real <lower=0> uAB;
  real <lower=0> uASC0;
  real <lower=0> sigma;
  real <lower=0> sigmapAB;
  real <lower=0> sigmauASC;
  //real <lower=0> rpABmu;
  //real <lower=0> ruASCmu;
  real <lower=0> rpAB[indivs];
  real <lower=0> ruASC[indivs];
}
model {
  vector[2] yhat[nobs];
  //prior distributions
  pAB0 ~ normal(0.5, 0.5);
  uASC0 ~ normal(0.5, 0.5);
  //rpABmu ~ normal(0, 0.001);
  uAB ~ normal(0, 0.001);
  //ruASCmu ~ normal(0, 0.001);
  sigmapAB ~ inv_gamma(0.01, 0.01);
  sigmauASC ~ inv_gamma(0.01, 0.01);
  sigma ~ inv_gamma(0.01, 0.01);
  for (j in 1:indivs) {
    real pAB = pAB0*rpAB[j];
    real uABt = uAB;
    real uASC = uASC0*ruASC[j];
    
    yhat = ode_rk45(model3, y0, t0, ts, pAB, uABt, uASC);
  //likelihood
  for (i in 1:nobs) {
    antib[i] ~ lognormal(yhat[i,2], sigma); 
    rpAB[j] ~ normal(0.01, sigmapAB);
    ruASC[j] ~ normal(0.01, sigmauASC);
  }
}
}

generated quantities {
  real z_pred[nobs];
  for (t in 1:nobs){
    z_pred[t] = lognormal_rng(antib[t], sigma);
  }
}
2

Hey,

no ODE expert here and I just have a quick guess: The location parameter of the Lognormal is on the log-scale, so maybe you need to do antib[i] ~ lognormal( log( yhat[i,2] ), sigma);? Right now Stan could be struggling with fairly large values for yhat (I suppose yhat is on the natural scale).

Cheers,
Max

3

Hi,

Thank you for your answer. I have tried this, but I got a even slower time:
Chain 1: Gradient evaluation took 0.03962 seconds
Chain 1: 1000 transitions using 10 leapfrog steps per transition would take 396.2 seconds.

4

In my opinion, the model does not look simple (there is a double for-loop inside with some sampling involved), so I may suspect slow convergence. For reference, sometimes I also have similar times for convergence and I think it may be normal for difficult problems

5

Hi,

When running the model I get the following warning message when I use cmdstanr and sample method:

warning
warning1042×469 36.8 KB

Could this be linked to slow convergence?

6

You can specify the initial values for the chains in order to avoid those warning messages

7

Hi, do you mean inside the sample function? Do you have any good references for this?

8

I definitely would not make a decision about normal versus lognormal based on the gradient evaluation time. This should be based on whether you expect the quantity on the left-hand-side to be normally or lognormally distributed.

My overall suggestion to you is to run a couple of chains in parallel and leave them alone for a couple of days until they finish. The output will be very informative about whether you are fitting the correct model (in which case the question will be whether there are any tweaks to be made to improve efficiency) versus fitting a poorly specified model or a model other than the one you intend to be fitting (in which case the question will be how to write down the model that you want). An ODE solve inside a for-loop is probably expected to take a while. This for-loop needs to execute hundreds of times per iteration, so you are dealing with tens of thousands of ODE solves per iteration.

9

All right, I’ll do that. Thanks a lot!

1 Like