Hi, I am triying to run this code in stan
smodel = "functions {
real[] sir(real t, real[] y, real[] theta,
real[] x_r, int[] x_i) {
real S = y[1];
real E = y[2];
real L = y[3];
real L1 = y[4];
real D = y[5];
real C2 = y[6];
real D2 = y[7];
real N = x_i[1];
// fixed parameters
real M = 4590490.56;
real mu = 1/(75.50);
real alphae = 0.00002;
real sigma1 = 0.01;
real epsilon = 1.02;
real rc = 0.00001;
real thetaest= 0.75;
real phi1= 0.05;
real pc;
if(t<18){
pc = 0;
}else{
pc = 0.65;
}
// parameters to be estimated
real nu =theta[1];
real tau = theta[2];
real sigma2 = theta[3];
real omega = theta[4];
real d1t= theta[5];
real betat= theta[6];
real lambdat = (betat*D)/N;
real dS_dt = M + alphae*E + sigma1*L + pc*sigma2*L1 + (thetaest*tau)*D - (mu + lambdat)*S;
real dE_dt = lambdat*(S + epsilon*L) - (mu+ nu + alphae)*E;
real dL_dt = (1- phi1)*nu*E +(1 - thetaest)*tau*D+(1-pc)*sigma2*L1 - (mu + sigma1 + rc + omega + epsilon*lambdat)*L;
real dL1_dt = rc*L - (mu + sigma2)*L1;
real dD_dt = phi1*nu*E + omega*L - (mu + d1t + tau)*D;
real dC2_dt = phi1*nu*E + omega*L;
real dD2_dt = d1t*D;
return {dS_dt, dE_dt,dL_dt,dL1_dt, dD_dt,dC2_dt, dD2_dt};
}
}
data {
int<lower=1> n_days;
real y0[7];
int to;
real t0;
real ts[n_days];
int N;
int cases1[to];
int cases2[to];
}
transformed data {
real x_r[0];
int x_i[1] = { N };
}
parameters {
// parameters to be estimated
real <lower =0> tau;
real <lower =0> sigma2;
real <lower =0> omega;
real <lower =0> nu;
real <lower =0> d1t;
real <lower =0> betat;
real<lower=0> phi_inv;
}
transformed parameters{
real theta[6];
real y[n_days, 7];
real incidence[n_days];
real deaths[n_days];
real phi = 1. / phi_inv;
{
theta[1] = nu;
theta[2] = tau;
theta[3] = sigma2;
theta[4] = omega;
theta[5] = d1t;
theta[6] = betat;
y = integrate_ode_rk45(sir, y0, t0, ts,theta, x_r, x_i);
incidence[1] = 92670;
deaths[1] = 5879;
for (i in 2:n_days){
incidence[i] = y[i, 6] - y[i-1, 6];
deaths[i] = y[i, 7] - y[i-1, 7];
}
}
}
model {
nu ~ normal (3,1) T[0,];
tau ~normal (0.80,0.15) T[0,];
sigma2 ~normal (1.2,0.10) T[0,];
omega ~normal (0.005,0.0015) T[0,];
d1t~normal (0.045,0.01) T[0,];
betat~normal (3.2,2) T[0,10];
phi_inv ~ exponential(5);
cases1 ~ neg_binomial_2(incidence[1:to], phi);
cases2 ~ neg_binomial_2(deaths[1:to], phi);
}
generated quantities {
real pred_cases1[n_days];
real pred_cases2[n_days];
// pred_cases1 = neg_binomial_2_rng(incidence, phi);
pred_cases1 = incidence;
//pred_cases2 = neg_binomial_2_rng(deaths, phi);
pred_cases2 = deaths;
}
"
The model runs well when I do not condition the parameter pc over time, but when I include this part (as the code above):
if(t<18){
pc = 0;
}else{
pc = 0.65;
}
the code throws me an error like this:
SYNTAX ERROR, MESSAGE(S) FROM PARSER:
Variable "real" does not exist.
31:
32: // parameters to be estimated
33: real nu =theta[1];
^
34: real tau = theta[2];
-------------------------------------------------
it does not look like a problem with the ordinary differential equations solver in STAN, Would you help me please?
Thank you in advance.