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 ?

``````data {
int nb_seq;
int n_pat;
real doseref;
vector[9] prior;
}
transformed data {
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){
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];
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)]);
}
}
}
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){
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.