Model estimation very long: how to gain time / optimize?

Dear all,

I am working on estimating the probability of a binary endpoint that we measure at different timepoints. Once it is observed at one timepoint, it is not longer followed and measured.
I build a stan that I fit with rstan, but it is very very long to run. I tried to optimize some parts, but I need to do better because I need to call this estimation a lot.
I would like to have your opinion and your help on how I could gain time when fitting it ?

Thank you in advance

data {
  int nb_seq;
  int nb_adm;
  int n_pat;
  int nb_adm_pat[n_pat];
  matrix[nb_seq,nb_adm] sequences;
  matrix[n_pat, nb_adm] seq_donnees;
  int y[n_pat,nb_adm];
  real doseref;
  vector[9] prior;
}
transformed data {
  matrix[nb_seq,nb_adm] norm_cum_sequences;
  matrix[nb_seq,nb_adm] norm_sequences;
  for(k in 1:nb_seq){
    norm_cum_sequences[k,] = cumulative_sum(sequences[k,]/doseref);
    norm_sequences[k,] = sequences[k,]/doseref;
  }
}
parameters {
  vector[3] param;
}
transformed parameters { 
  vector[3] exp_param;
  vector[nb_seq] ptox_seq;
  exp_param = exp(param);
  for(k in 1:nb_seq){
    vector[nb_adm] ptox_j;
    vector[nb_adm] ptox_cum_j;
    vector[nb_adm] temp_cumprod;
    real exp_calcul;
    exp_calcul = exp(param[1] - exp_param[2] + exp_param[3]*log(norm_sequences[k,1]));
    ptox_j[1] = exp_calcul / (1 + exp_calcul);
    ptox_cum_j[1] = ptox_j[1];
    temp_cumprod[1] = 1-ptox_j[1];
    for(j in 2:nb_adm){  
      exp_calcul = exp(param[1] - exp_param[2]*log(norm_cum_sequences[k,j-1]) + exp_param[3]*log(norm_sequences[k,j]));
      ptox_j[j] = exp_calcul / (1 + exp_calcul);
      temp_cumprod[j] = temp_cumprod[j-1] * (1-ptox_j[j]);
      ptox_cum_j[j] = ptox_j[1] + sum(ptox_j[2:j] .* temp_cumprod[1:(j-1)]);
    }
    ptox_seq[k] = ptox_cum_j[nb_adm];
  }
}
model {
  vector[3] mu;
  matrix[3,3] Sigma;
  mu[1] = prior[1];
  mu[2] = prior[2];
  mu[3] = prior[3];
  Sigma[1,1] = prior[4]*prior[4];
  Sigma[2,2] = prior[5]*prior[5];
  Sigma[3,3] = prior[6]*prior[6];
  Sigma[1,2] = prior[7]*prior[4]*prior[5];
  Sigma[2,1] = Sigma[1,2];
  Sigma[1,3] = prior[8]*prior[4]*prior[6];
  Sigma[3,1] = Sigma[1,3];
  Sigma[2,3] = prior[9]*prior[5]*prior[6];
  Sigma[3,2] = Sigma[2,3];
  param ~ multi_normal(mu, Sigma);
  for(i in 1:n_pat){
    for(j in 1:nb_adm_pat[i]){
      real ptox;
      real exp_calcul;
      if(j==1){
        exp_calcul = exp(param[1] - exp_param[2] + exp_param[3]*log(seq_donnees[i,1]/doseref));
      }
      else if(j>1){
        exp_calcul = exp(param[1] - exp_param[2]*log(sum(seq_donnees[i,1:(j-1)])/doseref) + exp_param[3]*log(seq_donnees[i,j]/doseref));
      }
      ptox = exp_calcul / (1 + exp_calcul);
      y[i,j] ~ bernoulli(ptox);
    }
  }
}

Hm. The way you’re passing in constants for use as “priors” is unusual; I suspect that such point-priors are uncommon for good reason.

Also, see here.