Hi all,
The following model takes 2.625 seconds, and it is taking really long to run.
I do not know if the problem belongs to the data.? (I have around 700 ids and 2085 obs), or the problem is the prior distributions? But I also do not think that the priors are very informative.
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];
real ts[N];
real t0;
}
transformed data {
real x_r[0];
int x_i[0];
}
parameters {
real logps;
real logpl;
real logmuab;
real logAB0;
real logps0;
real logpl0;
real <lower= 0> sigma;
}
model {
real new_mu[N,3];
real mu[1,3];
real eval_time[1];
real theta[3];
real inits[3];
real ps = exp(logps);
real pl = exp(logpl);
real muab = exp(logmuab);
real AB0 = exp(logAB0);
real ps0 = exp(logps0);
real pl0 = exp(logpl0);
//prior distributions
logps ~ normal(0, 1);
logpl ~ normal(0, 1);
logmuab ~ normal(0, 10);
logAB0 ~ normal(5, 10);
logps0 ~ normal(0, 1);
logpl0 ~ normal(0, 1);
sigma ~ normal(1, 10);
for (obs_id in 1:N){
inits[1] = AB0;
inits[2] = inits[1]/ps;
inits[3] = inits[1]/pl;
theta[1] = ps;
theta[2] = pl;
theta[3] = muab;
eval_time[1] = ts[obs_id];
if (eval_time[1] == 0){
new_mu[obs_id,1] = log(inits[1]) - pow(sigma, 2.0)/2;
new_mu[obs_id,2] = inits[2];
new_mu[obs_id,3] = inits[3];
}
else{
mu = integrate_ode_rk45(sldecay, inits, t0, eval_time, theta, x_r, x_i);
new_mu[obs_id,1] = log(mu[1,1]) - pow(sigma, 2.0)/2;
new_mu[obs_id,2] = mu[1,2];
new_mu[obs_id,3] = mu[1,3];
}
//likelihood function for AB levels
y[obs_id] ~ lognormal(new_mu[obs_id,1], 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);
}
}