Updating ode to new syntax. Exception: ode_rk45: ode paramers and data is inf

Modeling
, , ,
1

I’m trying to update some ode Stan code. My updated Stan code looks like this

functions {
  vector sir(real t, vector y, real beta, real gamma, array[] int x_i) {
    real S = y[1];
    real I = y[2];
    real R = y[3];
    real N = x_i[1];
      
    vector[3] dydt;	

    dydt[1] = -beta * I * S / N; // dS_dt
    dydt[2] =  beta * I * S / N - gamma * I; // dI_dt
    dydt[3] =  gamma * I; // dR_dt
      
    return dydt;
  }
}

data {
  int<lower=1> n_days;
  array[n_days] real t;
  vector[3] y0;
  int N;
  array[n_days] int cases;
}

transformed data {
  real t0 = 0;
  array[0] real x_r;
  array[1] int x_i = { N };
}

parameters {
  real<lower=0> beta;
  real<lower=0> gamma;
  real<lower=0> phi_inv_sqrt;
}

transformed parameters {
  real phi = 1. / (square(phi_inv_sqrt));	
  array[n_days] vector[3] y = ode_rk45(sir, y0, t0, t,
                                       beta, gamma, x_i);
}

model {
  beta ~ normal(1.7, 1.4); 
  gamma ~ normal(0.2, 0.53);
  phi_inv_sqrt ~ normal(0, 1);
	
  cases ~ neg_binomial_2(y[ : , 3], phi);
}

generated quantities {
  real R0 = beta / gamma;
  real recovery_time = 1 / gamma;
  array[n_days] real pred_cases;
  pred_cases = neg_binomial_2_rng(y[ : , 3], phi);
}

This runs but when I try to fit it I get many of these errors (from cmdstanr):
Chain 2 Informational Message: The current Metropolis proposal is about to be rejected because of the following issue:
Chain 2 Exception: ode_rk45: ode parameters and data is inf, but must be finite! (in ‘C:/…/.stan’, line 40, column 2 to line 41, column 57)
Chain 2 If this warning occurs sporadically, such as for highly constrained variable types like covariance matrices, then the sampler is fine,
Chain 2 but if this warning occurs often then your model may be either severely ill-conditioned or misspecified.

In total I get 185 of 2000 (9.0%) transitions ended with a divergence.

The relevant R code with data is given below.

library(outbreaks) # for dataset
library(rstan)

cases <- influenza_england_1978_school$in_bed 

# total count
N <- 763;

# times
n_days <- length(cases) 
t <- seq(0, n_days, by = 1)
t <- t[-1]

# initial conditions
i0 <- 1
s0 <- N - i0
r0 <- 0
y0 = c(S = s0, I = i0, R = r0)

# data for Stan
data_sir <- list(n_days = n_days, y0 = y0, t = t, N = N, 
                 cases = cases)

fit_sir_negbin <- sampling(model,
                           data = data_sir,
                           seed = 0)

I’m not sure where I messed up when trying to convert syntax of Stan model? I’m using rstan version 2.32.6 / cmdstanr version 2.35.0. I have attached the marginal posteriors I get (compare to posterior at the bottom from old code).

new_code_posterior
new_code_posterior861×717 80.7 KB

Old stan code
If its any help here is stan code using old syntax which runs fine.

// Code for ODE
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 dI_dt =  beta * I * S / N - gamma * I;
      real dR_dt =  gamma * I;
      
      return {dS_dt, dI_dt, dR_dt};
  }
}

data {
  int<lower=1> n_days;
  real y0[3];
  real t[n_days];
  int N;
  int cases[n_days];
}

// Track status of individuals in SIR model in x
transformed data {
  real t0 = 0; 
  real x_r[0];
  int x_i[1] = { N };
}

// Truncated prior using <lower=...>
parameters {
  real<lower=0> beta;
  real<lower=0> gamma;
  real<lower=0> phi_inv_sqrt;
}

transformed parameters{
  real y[n_days, 3];
  real phi = 1. / (square(phi_inv_sqrt));
  {
    real theta[2];
    theta[1] = beta;
    theta[2] = gamma;

    y = integrate_ode_rk45(sir, y0, t0, t, theta, x_r, x_i);
  }
}

model {
  //priors 
  beta ~ normal(1.7, 1.4); 
  gamma ~ normal(0.2, 0.53);
  phi_inv_sqrt ~ normal(0, 1);
  
  // sampling distribution
  // col(matrix x, int n) - The n-th column of matrix x. 
  // Here the number of infected people 
  cases ~ neg_binomial_2(col(to_matrix(y), 2), phi);
}

// Quantities of interest
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);
}

And here is the marginal posteriors I get from this code

old_code_posterior
old_code_posterior832×592 42.7 KB

2

Without combing through this line by line, I noticed you may be pulling different columns of y for the new syntax versus the old syntax. Not sure if this was intentional.

Here’s a version where I just replaced arrays with the new syntax, and tweaked the arguments for the ode:

// Code for ODE
functions {
  vector sir(real t, vector y, array[] real theta, 
             array[] real x_r, array[] 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 dI_dt =  beta * I * S / N - gamma * I;
      real dR_dt =  gamma * I;
      
      return [dS_dt, dI_dt, dR_dt]';
  }
}

data {
  int<lower=1> n_days;
  vector[3] y0;
  array[n_days] real t;
  int N;
  array [n_days] int cases;
}

// Track status of individuals in SIR model in x
transformed data {
  real t0 = 0; 
  array[0] real x_r;
  array[1] int x_i = { N };
}

// Truncated prior using <lower=...>
parameters {
  real<lower=0> beta;
  real<lower=0> gamma;
  real<lower=0> phi_inv_sqrt;
}

transformed parameters{
  array[n_days] vector[3] y;
  real phi = 1. / (square(phi_inv_sqrt));
  {
    array[2] real theta;
    theta[1] = beta;
    theta[2] = gamma;

    y = ode_rk45(sir, y0, t0, t, theta, x_r, x_i);
  }
}

model {
  //priors 
  beta ~ normal(1.7, 1.4); 
  gamma ~ normal(0.2, 0.53);
  phi_inv_sqrt ~ normal(0, 1);
  
  // sampling distribution
  // col(matrix x, int n) - The n-th column of matrix x. 
  // Here the number of infected people 
  cases ~ neg_binomial_2( y[,2] , phi);
}

// Quantities of interest
generated quantities {
  real R0 = beta / gamma;
  real recovery_time = 1 / gamma;
  array[n_days] real pred_cases;
  pred_cases = neg_binomial_2_rng( y[,2], phi);
}

It compiles fine, but I didn’t test it.

1 Like