ODE Errors & Problems with Initialization

Modeling
, , ,
1

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.

2

Branching on y_hat is not a good idea because the sampler expects the probability to be a continuous function of parameters. It’s probably better to add 1e-6 to all y_hat values unconditionally. And note the minus; it should not be 1e6 which is a very large number. This sign error is also present in y0 calculation.

1 Like
3

@nhuure: Thank you for catching that absolutely abysmal gaffe! That was a tremendous oversight.

I’ve removed the branching and it seems to be giving me my desired results. That I think was the key source of my trouble. Thank you for your help!

4

correct.

The ode integrator is only positive up to what is the absolute error tolerance. So just always add something which is bigger than that to the solution (always). Then it should be just fine as this is commonly only a tiny number compared to where your meaningful values live.

2 Likes