I am trying to pick up PyStan for my research project involving ODE. I have found this very instructive paper:

https://mc-stan.org/users/documentation/case-studies/boarding_school_case_study.html#2_using_simulated_data_to_understand_our_model

but they used Stan via R, so I am trying to reproduce part of their work using PyStan. Specifically my current goal is just to produce some prior predictive samples like the 3rd figure of Section 2 in the paper. I attach below the parameters I fed to the model and the âstanâ code that Iâve slightly modified to achieve my goal.

```
N = 763
t = np.array(range(1,31))
t0 = 0
i0 = 1
r0 = 0
y0 = [N-i0, i0, r0]
n_days = 30
simu_3 = dict(n_days = n_days, y0=y0, t0=t0, ts=t, N=N)
```

```
sir_stan = '''
functions {
real[] sir(real t, real[] y, real[] theta,
real[] x_r, int[] x_i) {
real S = y[1];
real I = y[2];
real R = y[3];
real N = x_i[1];
real beta = theta[1];
real gamma = theta[2];
real dS_dt = -beta * I * S / N;
real dR_dt = gamma * I;
real dI_dt = -(dS_dt + dR_dt);
return {dS_dt, dI_dt, dR_dt};
}
}
data {
int<lower=1> n_days;
real<lower=0> y0[3];
real t0;
real ts[n_days];
int N;
}
transformed data {
real x_r[0];
int x_i[1] = { N };
}
parameters {
real<lower=0> gamma;
real<lower=0> beta;
real<lower=0> phi_inv;
}
transformed parameters{
real y[n_days, 3];
real phi = 1. / phi_inv;
{
real theta[2];
theta[1] = beta;
theta[2] = gamma;
// Solve ODE at current parameter values
y = integrate_ode_rk45(sir, y0, t0, ts, theta, x_r, x_i);
}
// Adding a small +ve value to ensure strict positivity of z
for (k in 1:3){
for (n in 1:n_days){
y[n,k] = y[n,k] + 0.00001;
}
}
}
model {
//priors
beta ~ normal(2, 1);
gamma ~ normal(0.4, 0.5);
phi_inv ~ exponential(5);
}
generated quantities {
real R0 = beta / gamma;
real recovery_time = 1 / gamma;
real pred_cases[n_days];
pred_cases = neg_binomial_2_rng(col(to_matrix(y), 2), phi);
// col(to_matrix(y),2) means second column of the matrix y, that is I(t)
}
'''
```

Notice that Iâve added the part in the `transformed parameters`

block because I had some errors indicating that my input to the `negative_binomial_2_rng`

has negative location parameter, so I forced it to be positive in order to debug the problem. However when the sampling is done, and I proceed to check one of my sample of `y`

, I was expecting time series of S , I , R for 30 days. But the result was that, the values somewhat indicate that they âremained the sameâ throughout the 30 days and so I thought it was no wonder that it complained about negative values since âIâ never increased but stayed the same throughout the 30 days for all the samples.

Is there anything that I should further modify but that I missed to achieve my goal of getting prior predictive samples from this model? Thanks.