Dear all,
I have been stuck on this topic for some time now. I am trying to estimate time-varying infection (with other time-independent epidemiological) parameters for a complex multi-patch SEIRS model. But I decided to first implement the idea for a simple single-patch SEIRS model. I am not quite sure how to achieve this, but I imagine that I need to use a for loop at the transformed parameters block, to evaluate the ODE integrator at each time point. I tried implementing this, but on compiling, I got the following error:
86:
-------------------------------------------------
123: //params[4] = nu;
124:
125: y[1,] = ode_rk45(SEIRS, y0, -0.000001, dt[1],beta, kappa, gama, Lambda, mu, tau, phi); //integrate_ode_bdf when stiff
^
126:
127:
-------------------------------------------------
Ill-typed arguments supplied to function 'ode_rk45':
(<F1>, vector, real, real, vector, real, real, real, real, real, real)
where F1 = (real, vector, real, real, real, real, real, real, real) => vector
Available signatures:
(<F1>, vector, real, array[] real, real, real, real, real, real, real, real) => vector
The fourth argument must be array[] real but got real
I am using the new ODE interface, which requires that the times be given as array real. How to provide a single real time at each iteration is not very clear for me. I will appreciate it if you could help me debug the following Stan code, which I use cmstanpy to run.
functions {
array[] real diff(array[] real x){
int N = num_elements(x);
array[N-1] real result;
for (i in 1:(N-1)){
result[i] = x[i+1] - x[i];
}
return result;
}
vector SEIRS( real t, // time
vector y, // State S, E, I, R, N
real beta, // infection parameter beta
real kappa,
real gama,
real Lambda,
real mu,
real tau,
real phi
//vector n_obs
//vector n_fake
//real[] x_r, // real ODE argument that depend only on data
//int[] x_i
){ // integer ODE argument that depend only on data
vector[6] dydt;
dydt[1] = Lambda*y[6] - beta*y[1]*y[3]/y[6] - mu*y[1] + tau*y[4];
dydt[2] = beta*y[1]*y[3]/y[6] - (kappa + mu) * y[2];
dydt[3] = kappa * y[2] - (gama + phi + mu) * y[3];
dydt[4] = gama * y[3] - (tau + mu) * y[4];
dydt[5] = kappa * y[2];
dydt[6] = Lambda*y[6] - (mu * y[6] + phi * y[3]);
return dydt;
}
}
data{
int<lower = 1> n_obs; // Number of days sampled
int<lower = 1> n_fake; // number of fake days sampled for prediction
array[n_obs] int<lower = 0> cases; // The observed data (daily incidence)
vector<lower = 0>[6] y0; //initial conditions
real<lower = 0> t0; // Initial time point (zero)
array[n_obs+1] real<lower = t0> dt; // Time points that were sampled
array[n_fake] real<lower = t0> fake_ts; // Fake prediction time points that were sampled
array[n_obs] real<lower = 0> mob_data; //mobility data
real<lower = 0> N0; // Initial pop size
real<lower = 0> R0; //
real<lower = 0> Lambda; // We considere known
real<lower = 0> mu; // We considere known
real<lower = 0> tau; // We considere known
real<lower = 0> phi; // We considere known
}
//transformed data {
// real x_r[4] = {Lambda, mu, tau, phi};
// int x_i[0];
//}
parameters{
//Support of parameters
//real<lower = 0.5, upper = 2.0> beta_t;
vector<lower = 0.8, upper = 1.60>[n_obs] beta;
//vector<lower = 0.5, upper=1.6>[n_obs] beta0;
//real<lower = 0.5, upper=1.6> beta0;
//real<lower = 0.5, upper=1.6> beta1;
real<lower = 0.5, upper = 1.5> kappa;
real<lower = 0.5, upper = 1.5> gama;
//real<lower = 0.1> nu_inv;
real<lower = 0.1> nu_inv;
//real<lower = 0, upper = 15> E0;
//real<lower = 0, upper = 30> I0;
//real<lower = 5e-7> nu;
}
transformed parameters{
array[n_obs+1] vector<lower = 0.>[6] y; // Output from the ODE solver
array[n_fake] vector<lower = 0.>[6] y_pred_Inc; // Output from the ODE solver for prediction
real nu = 1./nu_inv;
//vector[n_obs+1] beta =
//{
// vector[6] y0; // initial conditions
// y0[1] = N0 - ( E0 + I0 + R0);
// y0[2] = E0;
// y0[3] = I0;
// y0[4] = R0;
//y0[5] = E0 + I0 + R0;
//y0[6] = N0;
{
//real params[3];
//params[1] = beta;
//params[2] = kappa;
//params[3] = gama;
//params[4] = nu;
y[1,] = ode_rk45(SEIRS, y0, -0.000001, dt[1],beta, kappa, gama, Lambda, mu, tau, phi); //integrate_ode_bdf when stiff
y_pred_Inc[1,] = ode_rk45(SEIRS, y0, -0.000001, fake_ts[1], beta, kappa, gama, Lambda, mu, tau, phi); //integrate_ode_bdf when stiff (for predictions)
for (t in 2:n_fake){
//real beta_t = exp(beta0[t] + beta1[t]*mob_data[t]);
real beta_t = beta[t]; //+ beta1*mob_data[t];
//real dt[1] = t - 1;
//real fake_ts[1] = t - 1;
//}
y[t,] = ode_rk45(SEIRS, y[t-1,], -0.000001, dt[t],beta_t, kappa, gama, Lambda, mu, tau, phi); //integrate_ode_bdf when stiff
// y0 = y[t];
y_pred_Inc[t,] = ode_rk45(SEIRS, y_pred_Inc[t-1,], -0.000001, fake_ts[t], beta_t, kappa, gama, Lambda, mu, tau, phi); //integrate_ode_bdf when stiff (for predictions)
// y0 = y_pred_Inc[t];
}
// }
}
}
model{
real epsilon = 1e-4;
// Prior distributions
for (t in 1:n_obs){
beta[t] ~ normal(-0.1115975, 0.2966405);
}
//beta0 ~ lognormal(-0.1115975, 0.2966405);
//beta1 ~ lognormal(-0.1115975, 0.2966405);
kappa ~ gamma(4.25173, 6.09586);
gama ~ gamma(4.25173, 6.09586);
//nu_inv ~ gamma(9.0045, 2.2508);
nu ~ gamma(0.64906, 0.28971);
//E0 ~ uniform(3.1699, 11.8301);
//I0 ~ uniform(6.3397, 23.6603);
// Likelihood
cases ~ neg_binomial_2(diff(y[ ,5]), nu+epsilon);
}
generated quantities {
real epsilon = 1e-4;
array[n_fake-1] real pred_Inc_cases;
pred_Inc_cases = neg_binomial_2_rng(diff(y_pred_Inc[ ,5]), nu+epsilon);
}