# How to speed up my Stan code?

**URL:** <https://discourse.mc-stan.org/t/how-to-speed-up-my-stan-code/22788>\
**Category:** Modeling\
**Tags:** fitting-issues, rstan\
**Created:** [May 31, 2021, 10:14am UTC](https://discourse.mc-stan.org/t/how-to-speed-up-my-stan-code/22788 "2021-05-31T10:14:32Z")\
**Posts on this page:** 9\
**Page:** 1

<div class="post-metadata">

**Author:** ![Isa\_Nazar](https://yyz2.discourse-cdn.com/flex030/user_avatar/discourse.mc-stan.org/isa_nazar/32/9376_2.png) [@Isa\_Nazar](https://discourse.mc-stan.org/u/Isa_Nazar)\
**Post date:** [May 31, 2021, 10:14am UTC](https://discourse.mc-stan.org/t/how-to-speed-up-my-stan-code/22788/1 "2021-05-31T10:14:32Z")

</div>

Hi Stan community,  
I’m new to Stan. The following paragraph is my stan code for the hierarchical survival model.  
I used this code for large data (N=35000).  
To run, It took ~10 days to even for 2000 iteration/500 warmup.  
I really want to know how can I speed up my code. So would be greatly appreciated if there are any tips!

Thank you.

My stan code is as below (N=35000, n\_grid=350, province=31, num=12).

fit ← stan(“final model.stan”,data = da,iter=iter, warmup=warmup, cores =4,chains = 2, thin = 1, init = “random”,control=list(max\_treedepth=15))

```stan
data {
  int <lower=0> province;
  int <lower=0> n_grid;
  int <lower=0> N;
  int <lower=0> num;
  int <lower=0> number;
  vector <lower=0> [province] p;
  matrix [province,n_grid] kernel;
  matrix [N,num] predictor;
  matrix [N,number] respon;
}
parameters{
  vector [num] beta;
  real <lower=0> sigma0;
  real <lower=0> sigma;
  vector [n_grid] x;
}
transformed parameters{
  real alpha = 1 / sigma;
  vector [N] lambdaa;
  vector[province] a; // `a` saved to output
     lambdaa = exp ((-1 * predictor * beta) / sigma);
           { // not part of output
    vector[province] z;
for (k in 1 : province) {
     vector[n_grid] landa_k;
for (j in 1 : n_grid) {
     landa_k[j] = kernel[k, j] * exp(x[j]);
           }
     z[k] = sum(landa_k) * p[k];
           }
     a = z / sum(z); // assign `a`
           }
           }
model{
      matrix [N,province] f;
      matrix [N,province] s;
     target += normal_lpdf(x|0,1);
    target += normal_lpdf(beta | 0, sigma0);
    target += inv_gamma_lpdf(sigma0| 0.001, 0.001);
    target += normal_lpdf(log (alpha) | 0, 5);
for (i in 1:N) {
for (k in 1:province){
    f[i,k]= a[k]*((lambdaa [i]*alpha*respon[i,1]^(alpha-1))/((1+(lambdaa [i]*respon[i,1]^alpha))^2));
    s[i,k]= a[k]*(1- ((lambdaa [i]*respon[i,1]^alpha)/ (1+(lambdaa [i]*respon[i,1]^alpha))));
    target += (respon[i,2] * log (f[i,k])) + ((1 - respon[i,2])* log(s[i, k]));
    }
  }
}

```

---

<div class="post-metadata">

**Author:** ![Ara\_Winter](https://yyz2.discourse-cdn.com/flex030/user_avatar/discourse.mc-stan.org/ara_winter/32/4284_2.png) [@Ara\_Winter](https://discourse.mc-stan.org/u/Ara_Winter)\
**Post date:** [June 10, 2021, 6:44pm UTC](https://discourse.mc-stan.org/t/how-to-speed-up-my-stan-code/22788/2 "2021-06-10T18:44:25Z")

</div>

Afternoon!

Could you provide the version of rstan, R, and your OS? Can you also provide a snippet of your data or some simulated data.

thanks

---

<div class="post-metadata">

**Author:** ![stevebronder](https://yyz2.discourse-cdn.com/flex030/user_avatar/discourse.mc-stan.org/stevebronder/32/17242_2.png) [@stevebronder](https://discourse.mc-stan.org/u/stevebronder)\
**Post date:** [June 10, 2021, 8:46pm UTC](https://discourse.mc-stan.org/t/how-to-speed-up-my-stan-code/22788/3 "2021-06-10T20:46:36Z")

</div>

I think you can calculate `z` with a matrix multiply and elementwise multiplication on p

```stan
data {
  int<lower=0> province;
  int<lower=0> n_grid;
  int<lower=0> N;
  int<lower=0> num;
  int<lower=0> number;
  vector<lower=0> [province] p;
  matrix[province,n_grid] kernel;
  matrix[N,num] predictor;
  matrix[N,number] respon;
}
parameters {
  vector[num] beta;
  real<lower=0> sigma0;
  real<lower=0> sigma;
  vector[n_grid] x;
}
transformed parameters{
  real alpha = 1 / sigma;
  vector[N] lambda;
  vector[province] a; // `a` saved to output
  lambda = exp((-1 * predictor * beta) / sigma);
  { // not part of output
    vector[province] z;
    /**
    for these sizes
    // matrix[province,n_grid] kernel;
    // vector[n_grid] x;
    I think you can just do a mat mul since
    the result dimensions are by province 
    // result[province] = kernel[province, n_grid] * x[n_grid]
    for (k in 1:province) {
      vector[n_grid] lambda_k;
      for (j in 1:n_grid) {
        lambda_k[j] = kernel[k, j] * exp(x[j]);
      }
    // then this just becomes an elementwise multiply
      z[k] = sum(lambda_k) * p[k];
    }
    
    */
    z = (kernel * exp(x)) .* p[k];
   // is this supposed to be a simplex?
    a = z / sum(z); // assign `a`
  }
}

model {
  matrix[N,province] f;
  matrix[N,province] s;
  target += normal_lpdf(x| 0, 1);
  target += normal_lpdf(beta | 0, sigma0);
  target += inv_gamma_lpdf(sigma0| 0.001, 0.001);
  target += normal_lpdf(log(alpha) | 0, 5);
  for (i in 1:N) {
    for (k in 1:province) {
      f[i, k]= a[k] * ((lambda[i] * alpha * respon[i, 1]^(alpha - 1)) / ((1 + (lambda[i] * respon[i, 1]^alpha))^2));
      s[i, k]= a[k] * (1 - ((lambda[i] * respon[i, 1]^alpha) / (1 + (lambda[i] * respon[i, 1]^alpha))));
      target += (respon[i, 2] * log(f[i, k])) + ((1 - respon[i, 2]) * log(s[i, k]));
    }
  }
}

```

In the future I’d recommend cleaning up spacing and indentation so it’s a little simpler to follow in the future. Do you have a paper associated with the model?

---

<div class="post-metadata">

**Author:** ![stevebronder](https://yyz2.discourse-cdn.com/flex030/user_avatar/discourse.mc-stan.org/stevebronder/32/17242_2.png) [@stevebronder](https://discourse.mc-stan.org/u/stevebronder)\
**Post date:** [June 10, 2021, 9:41pm UTC](https://discourse.mc-stan.org/t/how-to-speed-up-my-stan-code/22788/4 "2021-06-10T21:41:16Z")

</div>

> [@stevebronder](#):
>
> `(1 + (lambda[i] * respon[i, 1]^alpha)`

you can also pull out and vectorize operations like the above that are done multiple times

---

<div class="post-metadata">

**Author:** ![Isa\_Nazar](https://yyz2.discourse-cdn.com/flex030/user_avatar/discourse.mc-stan.org/isa_nazar/32/9376_2.png) [@Isa\_Nazar](https://discourse.mc-stan.org/u/Isa_Nazar)\
**Post date:** [June 11, 2021, 12:07pm UTC](https://discourse.mc-stan.org/t/how-to-speed-up-my-stan-code/22788/5 "2021-06-11T12:07:11Z")

</div>

I really appreciated your information and help. I’m working with R 4.0.3 and the last version of RSTAN.  
Actually, I made a lot of changes to the previous code and the latest version of my code is as below. Please check and let me know if you have any suggestions for improving its speed. I will be very grateful.

fit \< - stan(“final model.stan”,data = da,iter=iter, warmup=warmup, cores =4,chains = 4, thin = 1, init = “random”)

```stan
// My Stan codes for hierarchical survival model for large dataset(N=35000, n_grid=354, province=31, num=12).

data {
  int <lower=0> province;
  int <lower=0> n_grid;
  int <lower=0> N;
  int <lower=0> num;
  int <lower=0> number;
  int county[N];
  vector <lower=0> [province] p;
  //row_vector[n_grid] kernel[province];
  matrix [province,n_grid] kernel;
  matrix [N,number] respon;
}
parameters{
  vector [num] beta;
  real <lower=0> sigma;
  vector [n_grid] x;
}
transformed parameters{
  vector [num] time_ratio;
  vector [N] lambdaa;
  real alpha = 1 / sigma;
  vector [n_grid] exp_x;
  vector[province] a; // `a` saved to output
     time_ratio = exp(beta);
  exp_x = exp(x);
     lambdaa = exp ((-1 * respon[,6:17] * beta) / sigma);
{ // `z` not part of output
    vector[province] z;
    z = (kernel * exp_x) .* p;
    a = z / sum(z);
   // vector[province] landa;
    //for (k in 1 : province) {
      // landa[k] = kernel[k] * exp_x * p[k];
   // }
   //a = landa / sum(landa); // assign `a`
  }
  }
model{
  vector [N] f;
  vector [N] s;
  real log_alpha;
  log_alpha = log (alpha);
    target += normal_lpdf(x|0,5);
    target += normal_lpdf(beta | 0, pow(10,5));
    target += normal_lpdf(log_alpha | 0, 5);
for (i in 1:N) {
    f[i]= a[county[i]]*((lambdaa [i]*alpha*respon[i,4]^(alpha-1))/((1+(lambdaa [i]*respon[i,4]^alpha))^2));
    s[i]= a[county[i]]*(1- ((lambdaa [i]*respon[i,4]^alpha)/ (1+(lambdaa [i]*respon[i,4]^alpha))));
    target += respon[i,5] * (log (f[i])-log(s[i])) + log(s[i]) - log_alpha;
  }
}

```

---

<div class="post-metadata">

**Author:** ![stevebronder](https://yyz2.discourse-cdn.com/flex030/user_avatar/discourse.mc-stan.org/stevebronder/32/17242_2.png) [@stevebronder](https://discourse.mc-stan.org/u/stevebronder)\
**Post date:** [June 11, 2021, 6:25pm UTC](https://discourse.mc-stan.org/t/how-to-speed-up-my-stan-code/22788/6 "2021-06-11T18:25:39Z")

</div>

if f and s are being logged you can propogate that log through the model to make the multiplications and exponentiations into additions and multiplication

---

<div class="post-metadata">

**Author:** ![Isa\_Nazar](https://yyz2.discourse-cdn.com/flex030/user_avatar/discourse.mc-stan.org/isa_nazar/32/9376_2.png) [@Isa\_Nazar](https://discourse.mc-stan.org/u/Isa_Nazar)\
**Post date:** [June 11, 2021, 11:05pm UTC](https://discourse.mc-stan.org/t/how-to-speed-up-my-stan-code/22788/7 "2021-06-11T23:05:15Z")

</div>

Thank you for your help and reply. Now, I modified my code as below according to your instructions. I would greatly appreciate it if you could let me know when you have checked my code.

```stan
data {
  int <lower=0> province;
  int <lower=0> n_grid;
  int <lower=0> N;
  int <lower=0> num;
  int <lower=0> number;
  int county[N];
  vector <lower=0> [province] p;
  //row_vector[n_grid] kernel[province];
  matrix [province,n_grid] kernel;
  matrix [N,number] respon;
}
parameters{
  vector [num] beta;
  real <lower=0> sigma;
  vector [n_grid] x;
}
transformed parameters{
  vector [num] time_ratio;
  vector [N] lambdaa;
  real alpha = 1 / sigma;
  vector [n_grid] exp_x;
  vector[province] a; // `a` saved to output
     time_ratio = exp(beta);
  exp_x = exp(x);
     lambdaa = exp ((-1 * respon[,6:17] * beta) / sigma);
{ // `z` not part of output
    vector[province] z;
    z = (kernel * exp_x) .* p;
    a = z / sum(z);
   // vector[province] landa;
    //for (k in 1 : province) {
      // landa[k] = kernel[k] * exp_x * p[k];
   // }
   //a = landa / sum(landa); // assign `a`
  }
  }
model{
  // vector [N] f;
    vector [N] log_f;
  //vector [N] s;
   vector [N] log_s;
  real log_alpha;
  log_alpha = log (alpha);
    target += normal_lpdf(x|0,5);
    target += normal_lpdf(beta | 0, pow(10,5));
    target += normal_lpdf(log_alpha | 0, 5);
for (i in 1:N) {
 // f[i]= a[county[i]]*((lambdaa [i]*alpha*respon[i,4]^(alpha-1))/((1+(lambdaa 
    [i]*respon[i,4]^alpha))^2));
     log_f[i]= log(a[county[i]])+log(lambdaa[i])+log_alpha+alpha*log(respon[i,4])- 
     log(respon[i,4])- 2 * log(1+(lambdaa[i]*respon[i,4]^alpha));

// s[i]= a[county[i]]*(1- ((lambdaa [i]*respon[i,4]^alpha)/ (1+(lambdaa [i]*respon[i,4]^alpha))));
     log_s[i]=log(a[county[i]])+log(1- ((lambdaa [i]*respon[i,4]^alpha)/ (1+(lambdaa 
     [i]*respon[i,4]^alpha))));
   // target += respon[i,5] * (log (f[i])-log(s[i])) + log(s[i]) - log_alpha;
  }
     target += respon[,5] .* (log_f - log_s) + log_s - log_alpha;
}

```

---

<div class="post-metadata">

**Author:** ![stevebronder](https://yyz2.discourse-cdn.com/flex030/user_avatar/discourse.mc-stan.org/stevebronder/32/17242_2.png) [@stevebronder](https://discourse.mc-stan.org/u/stevebronder)\
**Post date:** [June 13, 2021, 11:24am UTC](https://discourse.mc-stan.org/t/how-to-speed-up-my-stan-code/22788/8 "2021-06-13T11:24:10Z")

</div>

I can’t check your code because I don’t have data associated with your model. I recommend following the SBC approach and building code for your data generating process and using that while you optimize the model to make sure you are getting the same results

---

<div class="post-metadata">

**Author:** ![Isa\_Nazar](https://yyz2.discourse-cdn.com/flex030/user_avatar/discourse.mc-stan.org/isa_nazar/32/9376_2.png) [@Isa\_Nazar](https://discourse.mc-stan.org/u/Isa_Nazar)\
**Post date:** [June 14, 2021, 11:38am UTC](https://discourse.mc-stan.org/t/how-to-speed-up-my-stan-code/22788/9 "2021-06-14T11:38:27Z")

</div>

Ok. Thanks a lot for your information and help.