Hi, I’m working off of a version of the source code shared: https://jrmihalj.github.io/estimating-transmission-by-fitting-mechanistic-models-in-Stan/

```
functions {
real[] SI(real t,
real[] y,
real[] params,
real[] x_r,
int[] x_i) {
real dydt[3];
dydt[1] = - params[1] * y[1] * y[2];
dydt[2] = params[1] * y[1] * y[2] - params[2] * y[2];
dydt[3] = params[2] * y[2];
return dydt;
}
}
data {
int<lower = 1> n_obs; // Number of days sampled
int<lower = 1> n_pop; // Number of hosts sampled at each time point.
int y[n_obs]; // The binomially distributed data
real t0; // Initial time point (zero)
real ts[n_obs]; // Time points that were sampled
int<lower = 1> n_pred; // This is to generate "predicted"/"unsampled" data
real pred_ts[n_pred]; // Time points for "predicted"/"unsampled" data
}
transformed data {
real x_r[0];
int x_i[0];
int n_params = 2; // Number of model parameters: beta, gamma
int n_difeq = 3; // Number of differential equations in the system
}
parameters {
real<lower = 0> params[n_params]; // Model parameters
real<lower = 0, upper = 1> S0; // Initial fraction of hosts susceptible
}
transformed parameters{
real<lower = -1e8, upper = 1> y_hat[n_obs, n_difeq]; // Output from the ODE solver
real<lower = 0, upper = 1> y0[n_difeq]; // Initial conditions for both S and I
y0[1] = S0 + 1e6;
y0[2] = 1 - S0 - 1e6;
y0[3] = 0;
y_hat = integrate_ode_rk45(SI, y0, t0, ts, params, x_r, x_i);
for (i in 1:n_obs) {
if (y_hat[i,2]<0)
{
y_hat[i, 2] = y_hat[i, 2] + 1e6;
}
}
}
model {
params ~ normal(0, 2); //constrained to be positive
S0 ~ normal(0.5, 0.5); //constrained to be 0-1.
y ~ binomial(n_pop, y_hat[, 2]); //y_hat[,2] are the fractions infected from the ODE solver
}
generated quantities {
// Generate predicted data over the whole time series:
real pred_I[n_pred, n_difeq];
real R_0;
pred_I = integrate_ode_rk45(SI, y0, t0, pred_ts, params, x_r, x_i);
R_0 = params[1] / params[2]; // Basic reproduction number
}
```

The code entails the creation of a mock SIR model for epidemiological modeling.

I use the following python code for initialization:

```
N = 500
I0 = 1
R0 = 0
S0 = N - (I0 + R0)
y0 = (S0, I0, R0)
t = np.linspace(1, 700, 50)
gamma = 1.0/14.0
beta = 1.5 * gamma
pred_t = np.linspace(700, 1000, 15)
# to create dummy test data
y = odeint(func=sir_model, y0=y0, t=t, args=(N, beta, gamma))
data_2 = {
'n_obs' : len(t),
'ts' : t,
'y' : y,
'n_pop' : N,
't0' : 0.0,
'n_pred': len(pred_t),
'pred_ts' : pred_t
}
sm2 = pystan.StanModel(
model_name='sir_model_2',
model_code=stan_model_2)
fit_2 = sm2.sampling(data=data_2, iter=4000, warmup=1000, chains=2, n_jobs=-1, init='0')
```

I run into the following issue:

```
Rejecting initial value:
Error evaluating the log probability at the initial value.
Exception: validate transformed params: y_hat[i_0__][i_1__] is nan, but must be greater than or equal to -1e+08 (in 'unknown file name' at line 49)
RuntimeError: Initialization failed.
```

I earlier faced an issue of:

```
Rejecting initial value:
Error evaluating the log probability at the initial value.
Exception: binomial_lpmf: Probability parameter[5] is -4.31213e-07, but must be in the interval [0, 1] (in 'unknown file name' at line 63)
```

which prompted me to add the `1e6' and `

1e8` conditions.