So my latest update is as follow:

I am following the link shared by @bbbales2. My data contains daily number of confirmed cases of covid-19 in the UK. My stan code is as follow:

```
functions {
real[] sir(real t,
real[] y,
real[] theta,
real[] x_r,
int[] x_i) {
// I renamed my ODEs variables here:
real S = y[1];
real E = y[2];
real G = y[3];
real P = y[4];
real I = y[5];
real A = y[6];
real B = y[7];
real C = y[8];
real N = x_i[1];
// Here are the list of my parameters
real phi = theta[1];
real pe = theta[2];
real pie = theta[3];
real pee = theta[4];
real beta = theta[5];
real mu = theta[6];
real betam = theta[7];
real tau = theta[8];
real alpha = theta[9];
real Ts = theta[10];
real Ta = theta[11];
real Tm = theta[12];
real Ti = theta[13];
// here are the ODEs. I notice that in the link they define them as real dS_dt rather than dydt[1]
real dS_dt = - mu * beta * S * A - beta * S * I - beta * C * S - betam * B * S;
real dE_dt = mu * beta * S * A + (1 - phi) * beta * S * I + (1 - phi) * beta * C * S + (1 - phi) * betam * B * S - E / tau;
real dG_dt = pe * phi * beta * S * I + pie * phi * beta * C * S + pee * phi * betam * B * S - G / tau;
real dP_dt = (1 - pe) * phi * beta * S * I + (1 - pie) * phi * beta * C * S + (1 - pee) * phi * betam * B * S - P / tau;
real dI_dt = alpha * E / tau - I / Ts;
real dA_dt = (1 - alpha) * E / tau - A / Ta;
real dB_dt = G / tau - B / Tm;
real dC_dt = P / tau - C / Ti;
return {dS_dt, dE_dt, dG_dt, dP_dt, dI_dt, dA_dt, dB_dt, dC_dt};
}
}
data {
int<lower=1> n_days;
real y0[8];
real t0;
real ts[n_days];
int N;
int cases[n_days];
}
transformed data {
real x_r[0];
int x_i[1] = { N };
}
parameters {
real<lower = 0, upper = 1> phi;
real<lower = 0, upper = 1> pe;
real<lower = 0, upper = 1> pie;
real<lower = 0, upper = 1> pee;
real<lower = 0> beta;
real<lower = 0> mu;
real<lower = 0>betam;
real<lower = 0> tau;
real<lower = 0> alpha;
real<lower = 0> Ts;
real<lower = 0> Ta;
real<lower = 0> Tm;
real<lower = 0> Ti;
}
transformed parameters{
real y_hat[n_days, 8];
real<lower=0, upper=1> theta1[4];
real<lower = 0> theta2[9];
y_hat = integrate_ode_rk45(sir, y0, t0, ts, {phi, pe, pie, pee, beta, mu, betam, tau, alpha, Ts, Ta, Tm, Ti}, x_r, x_i);
for(i in 1:n_days){
for(j in 1:8){
y_hat[i, j] = 1e-6 + (1-2e-6)*y_hat[i, j];
}
}
}
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
cases ~ binomial(col(to_matrix(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));
}
```

But I can not run this again, as it says:

`Error in close.connection(zz) : cannot close 'message' sink connection`

when compiling the stan file. My R code is here along with other files. I wonder what am I doing wrong? any help is much appreciated. Could this be a problem at model block? I am defining the likelihood as binomial for infections which is the 5th ODE.

I also wonder why am I getting error as above only? while others seem to get pinpointed on where the error is happening?