Recently I have opened a topic here seeking a solution to the error I was getting while running stan. So far, I have managed to polish my code using the advice I got from that topic. Here is my stan code:

```
cat(
'
functions {
real[] SI(real t,
real[] y,
real[] params,
real[] x_r,
int[] x_i) {
real dydt[8];
real phi = params[1];
real pe = params[2];
real pie = params[3];
real pee = params[4];
real beta = params[5];
real mu = params[6];
real betam = params[7];
real tau = params[8];
real alpha = params[9];
real Ts = params[10];
real Ta = params[11];
real Tm = params[12];
real Ti = params[13];
dydt[1] = - mu * beta * y[1] * y[6] - beta * y[1] * y[5] - beta * y[8] * y[1] - betam * y[7] * y[1];
dydt[2] = mu * beta * y[1] * y[6] + (1 - phi) * beta * y[1] * y[5] + (1 - phi) * beta * y[8] * y[1] + (1 - phi) * betam * y[7] * y[1] - y[2] / tau;
dydt[3] = pe * phi * beta * y[1] * y[5] + pie * phi * beta * y[8] * y[1] + pee * phi * betam * y[7] * y[1] - y[3] / tau;
dydt[4] = (1 - pe) * phi * beta * y[1] * y[5] + (1 - pie) * phi * beta * y[8] * y[1] + (1 - pee) * phi * betam * y[7] * y[1] - y[4] / tau;
dydt[5] = alpha * y[2] / tau - y[5] / Ts;
dydt[6] = (1 - alpha) * y[2] / tau - y[6] / Ta;
dydt[7] = y[3] / tau - y[7] / Tm;
dydt[8] = y[4] / tau - y[8] / Ti;
return dydt;
}
}
data {
int<lower = 1> n_obs; // Number of days sampled
int<lower = 1> n_params; // Number of model parameters
int<lower = 1> n_difeq; // Number of differential equations in the system
int<lower = 1> n_sample;
int y[n_obs]; // The binomially distributed data
real t0; // Initial time point (zero)
real ts[n_obs]; // Time points that were sampled
}
transformed data {
real x_r[0];
int x_i[0];
}
parameters {
real<lower=0, upper=1> params1[4];
real<lower = 0> params2[9];
real<lower = 0, upper = 1> S0; // Initial fraction of susceptible
}
transformed parameters{
real y_hat[n_obs, n_difeq]; // Output from the ODE solver
real y0[n_difeq]; // Initial conditions for both S and I
y0[1] = 1 - 0.02 - 0.001;
y0[2] = 0.02;
y0[3] = 0;
y0[4] = 0;
y0[5] = 0.001;
y0[6] = 0;
y0[7] = 0;
y0[8] = 0;
y_hat = integrate_ode_rk45(SI, y0, t0, ts, append_array(params1, params2), x_r, x_i);
y_hat = 1e-6 + (1-2e-6)*y_hat;
}
model {
phi ~ normal(0, 1); //phi
pe ~ normal(0, 1); //pe
pie ~ normal(0, 1); //pie
pee ~ normal(0, 1); //pee
beta ~ normal(0.5, 0.9); //beta
mu ~ normal(0.2, 0.55); //mu
betam ~ normal(0, 1); //beta_m
tau ~ normal(1, 20); //tau
alpha ~ normal(0.1, 0.99); //alpha
Ts ~ normal(3, 20); //Ts
Ta ~ normal(1, 20); //Ta
Tm ~ normal(1, 20); //Tm
Ti ~ normal(0.5, 8); //Ti
S0 ~ normal(0.5, 0.5); //constrained to be 0-1.
y ~ binomial(n_sample, y_hat[,5]); //y_hat[,5] are the fractions infected from the ODE solver
}
generated quantities {
real R_0;
R_0 = beta * Ta * (alpha * Ts * (1 - phi) / Ta + mu * (1 - alpha));
}
',
file = "SI_fit.stan", sep="", fill=T)
```

As you can see I have 8 ODEs and 13 parameters to estimate, 4 of these parameters namely `phi`

, `pe`

, `pie`

and `pee`

are probabilities and so should be between 0 and 1. I was suggested to distinguish my parameters - those that are probabilities and those which are not, so I now have:

```
parameters {
real<lower=0, upper=1> params1[4];
real<lower = 0> params2[9];
...
}
```

and then in transformed parameters my y_hat is:

```
transformed parameters{
...
y_hat = integrate_ode_rk45(SI, y0, t0, ts, append_array(params1, params2), x_r, x_i);
y_hat = 1e-6 + (1-2e-6)*y_hat;
}
```

where I am appending params1 (probabilities) and params2 (the rest of the parameters). As you can see in the stan code now I have 2 names for my parameters: `params`

then in parameters I separate them into `params1`

and `params2`

I wonder if this makes any issue on running the stan code as it simply stucks at running `stan()`

. I am adding my whole code to this topic so one can run it at ease. Any help would be much appreciated by this stan newbie.

I am also dubious on my definition of model on R I defined:

```
stan_d = list(n_obs = N,
n_params = 13,
n_difeq = 8,
n_sample = 66650000,
y = cases,
t0 = 0,
ts = sample_time)
```

where N is the length of my list in days, I wonder if I need the n_sample which is the total population?

Issue1.R (4.6 KB)