Initialization Error in Stochastic Epidemic SEIR Model

[maxbiostat]

I was able to fit the model by replacing the binomial with a continuous approximation. I make no claims of correctness or appropriateness. Program is below.

functions{
  real cont_binomial_lpdf(real x, real n, real theta){
    real ans = lchoose(n, x)  + x*log(theta) + (n-x)*log1m(theta);
    if(x > n) return(negative_infinity());
    return(ans);
  }
}
data {
  int<lower=1> n_days;// total # of time point 
  int S0[n_days]; //of susceptible individuals at the start of the epidemic 
  int E0[n_days]; //of exposed individuals at the start of the epidemic
  int I0[n_days]; //of infectious individuals at the start of the epidemic
  int R0[n_days]; //of removed individuals at the start of the epidemic
  int N;          //population size
  int C[n_days];  // empirical data
  int D[n_days];  // empirical data
  int tcon;       // time point where control measures are introduced 
  int B0[n_days]; // initial conditions for unobserved entity. 
}
parameters {
  real<lower=0,upper=1> beta;
  real<lower=1/10.5,upper=1/3.5> gamma;
  real<lower=1/21.0,upper=1> rho;
  real<lower=0> q;
  real<lower=0> B[n_days]; // unobserved entity. # of susceptitble people who become infected
}
transformed parameters{
  real<lower=0> S[n_days]; // # of susceptible individuals
  real<lower=0> E[n_days]; // # of exposed individuals
  real<lower=0> I[n_days]; // # of infectious individuals
  real R[n_days]; // # of removed individuals
  real<lower=0> P[n_days]; //  probability of leaving S compartment
  real<lower=0> Beta[n_days];
  real<lower=0, upper=1> pc; //probability of leaving E compartment
  real<lower=0, upper=1> pr; //probability of leaving I compartment
  
  pc = 1 - exp(-rho);
  pr = 1- exp(-gamma);
  
  S[1] = S0[1] - B0[1];
  E[1] = E0[1] + B0[1] - C[1];
  I[1] = I0[1] + C[1] - D[1];
  R[1] = N - S[1] - E[1] - I[1];
  Beta[1] =  beta*exp(q*tcon); // dubious, please check
  P[1] = 1 - exp(-(Beta[1]*I[1]/N));
  
  for(t in 2:n_days){
    if(t >= tcon){
      Beta[t] = beta*exp(-q*(t-tcon));
    }
    else{
      Beta[t] = beta;
    }
    S[t]= S[t-1] - B[t-1];
    E[t]= E[t-1] + B[t-1] - C[t-1];
    I[t]= I[t-1] + C[t-1] - D[t-1];
    P[t] = 1 - exp(-(Beta[t]*I[t]/N));
    R[t] = N - S[t] - E[t] - I[t];
  }
  // print(R);
}
model{
  beta ~ gamma(2,10); // prior
  gamma ~  gamma(2,14); // prior
  rho ~  gamma(2,10); // prior
  q ~ gamma(2,10); // prior
  for(t in 2:n_days){
    B[t] ~ cont_binomial(S[t], P[t]); //  
    C[t] ~ cont_binomial(E[t], pc); // # of cases by date of symptom onset
    D[t] ~ cont_binomial(I[t], pr); // # of cases removed from infectious class
  }
}