Hi,

When I run the following ODEs with random effects, it takes too long, even if it says that it takes 0.013 seconds, which I do not see long.

Previously I ran this model without random effects, in which I called the ode function in the transformed parameters (without any loop), where I also specified the rest of the parameter transformations and was working well.

Now that I have the random effects and I have to define it in a loop (to specify them for each individual), I need to define this in the model statement, and therefore I had to move everything in the model statement.

What could improve the memory of this code?

```
functions {
real[] sldecay(real t, real[] y, real[] theta, real[] x_r, int[] x_i) {
real dydt[3];
dydt[1] = -theta[1]*y[1];
dydt[2] = -theta[2]*y[2];
dydt[3] = theta[1]*y[1] + theta[2]*y[2] - theta[3]*y[3];
return dydt;
}
}
data {
int<lower=1> N;
int<lower=1> nsubjects;
int<lower=1> subj_ids[N];
real <lower=0> y[N];
int <lower=1> ntimes;
real ts[ntimes];
real t0;
}
transformed data {
real x_r[0];
int x_i[0];
}
parameters {
real logps;
real logpl;
real logmuab;
real <lower = 0> sigmalogmuab;
real logAB0;
real <lower = 0> sigmalogAB0;
real logps0;
real logpl0;
real <lower= 0> sigma;
real rho_rhobit;
vector[2] r[nsubjects];
}
model {
real inits[3];
real mu[N, 3];
real theta[3];
real new_mu[N,3];
vector[nsubjects] rAB0;
vector[nsubjects] rmuab;
// transformed parameters
real ps = exp(logps);
real pl = exp(logpl);
real muab = exp(logmuab);
real sigma_muab = exp(sigmalogmuab);
real AB0 = exp(logAB0);
real sigma_AB0 = exp(sigmalogAB0);
real ps0 = exp(logps0);
real pl0 = exp(logpl0);
// defining the correlation
vector[2] mean_vec;
matrix[2,2] Sigma;
real rho = (exp(rho_rhobit) - 1)/(exp(rho_rhobit) + 1);
mean_vec[1] = -pow(sigma_AB0, 2.0)/2;
mean_vec[2] = -pow(sigma_muab, 2.0)/2;
Sigma[1,1] = pow(sigma_AB0, 2.0);
Sigma[2,2] = pow(sigma_muab, 2.0);
Sigma[1,2] = rho*sigma_AB0*sigma_muab;
Sigma[2,1] = rho*sigma_AB0*sigma_muab;
//prior distributions
logps ~ normal(0, 1);
logpl ~ normal(0, 1);
logmuab ~ normal(0, 10);
sigmalogmuab ~ normal(0, 1);
logAB0 ~ normal(5, 10);
sigmalogAB0 ~ normal(0, 1);
logps0 ~ normal(0, 1);
logpl0 ~ normal(0, 1);
sigma ~ normal(1, 10);
rho_rhobit ~ uniform(-10,10);
theta[1] = ps;
theta[2] = pl;
// defining the distribution for the random effects
for (subj_id in 1:nsubjects) {
r[subj_id] ~ multi_normal(mean_vec, Sigma);
rAB0[subj_id] = r[subj_id,1];
rmuab[subj_id] = r[subj_id,2];
}
for (obs_id in 1:N){
inits[1] = AB0*exp(rAB0[subj_ids[obs_id]]);
inits[2] = inits[1]/ps0;
inits[3] = inits[1]/pl0;
theta[3] = muab*exp(rmuab[subj_ids[obs_id]]);
}
//defining ode_solver
mu = integrate_ode_rk45(sldecay, inits, t0, ts, theta, x_r, x_i);
for (obs in 1:N){
new_mu[obs,1] = mu[obs,1];
new_mu[obs,2] = mu[obs,2];
new_mu[obs,3] = log(mu[obs,3]) - pow(sigma, 2.0)/2;
//likelihood function for AB levels
y[obs] ~ lognormal(new_mu[obs,3], sigma);
}
}
generated quantities {
real z_pred[N];
real sigma2 = sigma;
for (t in 1:N){
z_pred[t] = lognormal_rng(log(y[t] + 1e-5), sigma2);
}
}
```