Dear Stan-community,

I have an ODE with some fixed parameters. I actually have to parametrize three parameters (beta and kappa, p_reported, below code) from the fitting by using real infectious disease data. I am getting errors in the initial value although I defined the lower bound for them.

Rejecting initial value:

Error evaluating the log probability at the initial value.

Exception: neg_binomial_2_lpmf: Location parameter[…] is -40725.6, but must be > 0!

Any hint will be appreciated.

PS. Please don’t recommend using distribution for these fixed parameters. It is tried and will not solve the problem. Also, I prefer to keep the number of unknown parameters low.

```
functions {
real[] sir(real t, real[] y, real[] theta,
real[] x_r, int[] x_i) {
int N = x_i[1];
real beta = theta[1];
real kappa =theta[2];
real ia0 = theta[3];
real ip0 = theta[4];
real ic0 = theta[5];
real e0 = theta[6];
real a = x_r[1];
real ca = x_r[2];
real alpha= x_r[3];
real init[6] = {N - e0- ia0 -ip0 -ic0, e0, ia0, ip0, ic0, 0}; // initial values
real S = y[1] + init[1];
real E = y[2] + init[2];
real Ia = y[3] + init[3];
real Ip = y[4] + init[4];
real Ic = y[5] + init[5];
real R = y[6] + init[6];
real dS_dt = -beta * (ca*Ia + Ip + Ic) * S / N;
real dE_dt = beta * (ca*Ia + Ip + Ic) * S / N - a * E;
real dIa_dt = alpha* a * E - ca* Ia + kappa *Ia;
real dIp_dt = (1 -alpha)* a *E - Ip;
real dIc_dt = Ip - Ic;
real dR_dt = ca* Ia + Ic ;
return {dS_dt, dE_dt, dIa_dt, dIp_dt, dIc_dt, dR_dt};
}
}
data {
int<lower=1> n_days;
real t0;
real ts[n_days];
int N;
int cases[n_days];
}
transformed data {
int x_i[1] = {N};
real a = 0.30;
real ca = 0.66;
real alpha = 0.33;
real x_r[3];
x_r[1] = a;
x_r[2] = ca;
x_r[3] = alpha;
}
parameters {
real<lower=0> beta;
real<lower=0> kappa;
real<lower=0> phi_inv;
real<lower=0, upper=1> p_reported;
real<lower=0> ia0;
real<lower=0> ip0;
real<lower=0> ic0;
real<lower=0> e0;
}
transformed parameters{
real y[n_days, 6];
real incidence[n_days - 1];
real phi = 1. / phi_inv;
real theta[6];
theta = {beta,kappa,ia0,ip0,ic0,e0};
y = integrate_ode_rk45(sir, rep_array(0.0, 6), t0, ts, theta, x_r, x_i);
for (i in 1:n_days-1){
//incidence[i] = -(y[i+1, 4] - y[i, 4]+ y[i+1, 3] - y[i, 3]+ y[i+1, 2] - y[i, 2] + y[i+1, 1] - y[i, 1])* p_reported;
incidence[i] = (y[i+1, 5] - y[i, 5]) * p_reported;
}
}
model {
//priors
beta ~ normal(2, 0.5);
phi_inv ~ exponential(3);
p_reported ~ beta(1, 2);
ia0 ~ normal(2, 0.1);
ip0 ~ normal(2,0.1);
ic0 ~ normal(0.4,0.5);
e0 ~ normal(0.4,0.5);
kappa ~ normal(10, 0.1);
//sampling distribution
//col(matrix x, int n) - The n-th column of matrix x. Here the number of infected people
cases[1:(n_days-1)] ~ neg_binomial_2(incidence, phi);
}
generated quantities {
real pred_cases[n_days-1];
pred_cases = neg_binomial_2_rng(incidence, phi);
}
```