Multivariate ODE initialization issues

Modeling
, , ,
1

I am new to coding models directly in STAN and this also happens to be my first attempt at fitting an ODE model using it. I am having a few issues with the model initialization and have not been able to gather from the forums a solution that I can generalize to fix my issues. I receive the following errors when running the model for all chains and the sampling does not occur:

“Chain 1: Rejecting initial value:
Chain 1: Error evaluating the log probability at the initial value.
Chain 1: Exception: normal_lpdf: Location parameter is nan, but must be finite! (in ‘string’, line 80, column 2 to column 45)”

Using print() I have found none of the inputs to the ODE integrator are NA or unexpected but the matrix returned, ‘y_sim’ in the model, returns all “nan”.

I am running the model in rstan v2.32.5. Any help on debugging or fixing this error, or comments on the model in general would be much appreciated.

Thanks.

Relevant example below.

functions {
  real[] C_R(real t,           // Time
             real[] y,         // State variables: y[1] = R, y[2] = C
             real[] theta,     // Parameters
             real[] x_r,       // Real data (temperature, light)
             int[] x_i) {     // Int data (empty here)
    // observed states
    real R = y[1];
    real C = y[2];
    real cB = y[3];
    // parameters
    real a = theta[1];
    real h = theta[2];
    real r = theta[3];
    real K = theta[4];
    real e = theta[5];
    real m = theta[6];
    // real data

    real dcB_dt = a*R*C/(1+a*h*R);    // Biomass consumption
    real dR_dt = r*R*(1 - R/K) - cB;  // Prey dynamics
    real dC_dt = (e * cB) - (m * C);  // Predator dynamics
    return {dR_dt, dC_dt, dcB_dt};
  }
}

data {
  int<lower=1> T;             // Number of time steps
  real<lower=0> ts[T];        // Time values
  real<lower=0> R_initial;    // Initial resource population
  real<lower=0> C_initial;    // Initial consumer population
  real<lower=0> cB_initial;   // Initial consumption
  real<lower=0> C_obs[T];     // Observed mean consumer population data
  real<lower=0> R_obs[T];     // Observed mean resource population data
  real<lower=0> cB_obs[T];    // Observed Biomass specific consumption
  real<lower=0> sigma_C[T];   // Observed sd consumer population data
  real<lower=0> sigma_R[T];   // Observed sd resources population data
  real<lower=0> sigma_cB[T];  // Observed sd consumption
  // real<lower=0> temp[T];      // observed temperature data
  // real<lower=0> light[T];     // observed light data
}
transformed data{
  real x_r[0];
  int x_i[0];
}
parameters {
  real<lower=0,upper=1> a;         // attack/clearance rate parameter
  real<lower=0> h;         // handling time
  real<lower=0> r;         // prey growth rate
  real<lower=0> K;         // Carrying capacity
  real<lower=0,upper=1> m;         // Consumer mortality rate
  real<lower=0,upper=1> e; // gross growth efficiency
  // real<lower=0> sigma;
}

model {
  real y_init[3] = {R_initial, C_initial, cB_initial};
  real theta[6] = {a, h, r, K, e, m};  // Parameters
  // matrix[T, 2] y_sim;  // Simulated population dynamics
  real y_sim[T,3];  // Simulated population dynamics
  
  // Priors
  r ~ lognormal(1, 1);
  K ~ lognormal(10, 3);
  a ~ beta(1, 1);
  h ~ lognormal(0, 1);
  e ~ beta(1, 4);
  m ~ beta(1, 1);
  // sigma ~ exponential(0.5);
  print(a);
  print(K);
  print(r);
  print(h);
  print(e);
  print(m);
  
  // Solve ODE
  // y_sim = integrate_ode_bdf(C_R, y_init, 0, ts, theta, x_r, x_i);
  y_sim = integrate_ode_rk45(C_R, y_init, 0, ts, theta, x_r, x_i);
  print(y_sim);

// Likelihood
for (t in 1:T) {
  R_obs[t] ~ normal(y_sim[t, 1], sigma_R[t]);
  C_obs[t] ~ normal(y_sim[t, 2], sigma_C[t]);
  cB_obs[t] ~ normal(y_sim[t,3], sigma_cB[t]);
}

// for (t in 1:T) {
//   R_obs[t] ~ normal(y_sim[t, 1], sigma);
//   C_obs[t] ~ normal(y_sim[t, 2], sigma);
// }
}

generated quantities {
  real y_init[3] = {R_initial, C_initial, cB_initial};
  real theta[6];
  real y_pred[T, 3];
  
  // Parameters
  theta[1] = a;            // Saturation parameter
  theta[2] = h;            // Half-saturation constant (for resource)
  theta[3] = r;            // Base prey growth rate
  theta[4] = K;            // Carrying capacity
  theta[5] = e;            // Error term for consumer population
  theta[6] = m;            // Consumer mortality rate

  // Solve ODE for generated quantities
    // y_pred = integrate_ode_bdf(C_R, y_init, 0, ts, theta, x_r, x_i);
    y_pred = integrate_ode_rk45(C_R, y_init, 0, ts, theta, x_r, x_i);

}

The data I have been passing to this model are:

stan_data <- list(T = 10L, ts = c(41, 87, 142, 171, 208, 214, 259, 296, 317, 
346), R_initial = 7.95909938, C_initial = 317.44602496823, cB_initial = 2046.65868443038, 
    R_obs = c(7.95909938, 63.582298132, 166.26195652, 203.014871545303, 
    249.90652175, 269.91180124, 401.51878882, 146.994720504, 
    163.72037267184, 186.817701856), sigma_R = c(8.35977921339248, 
    32.7304043315385, 40.4895905956132, 77.1965146978821, 124.029486828363, 
    113.278703479259, 139.435933432936, 80.8086875268657, 99.999220067294, 
    126.500431670743), C_obs = c(317.44602496823, 925.57349071122, 
    5755.67593050651, 6409.7857725782, 15814.0416956033, 2492.23227708924, 
    2061.77271128308, 677.379879132129, 167.537343626406, 305.172916746582
    ), sigma_C = c(131.355718609556, 227.261418700627, 932.462821358444, 
    1316.3381507101, 6316.15110556263, 300.928146091724, 377.661169342442, 
    145.108450308571, 46.0375633029446, 71.6065909214821), cB_obs = c(2046.65868443038, 
    2060.03573956124, 376.377453805423, 763.560069928162, 33.6712991904767, 
    604.310412430487, 525.789893863879, 351.698498991322, 1655.07227107854, 
    2265.43247664591), sigma_cB = c(108.54143636971, 62.5486724534629, 
    7.96306161036773, 24.8668267776185, 0.933383933755735, 16.0339488862616, 
    11.4392829893471, 6.27473216475427, 60.7508978815252, 80.4613467879306
    ))
2

I had an issue like this somewhat recently, where under certain parameter combinations, the dynamics of the ODE exploded and quickly pushed some of the state variables to NaN.
I diagnosed this by placing a conditional reject() statement in the ODE function. You could add something like this to the end of your ode function:

if(is_nan(dR_dt) || is_nan(dC_dt) || is_nan(dcB_dt)) { 
    reject("NaN detected; t = ", t, "; derivs = ", [dR_dt, dC_dt, dcB_dt], 
              '; y = ", y, "; theta = ", theta);
return {dR_dt, dC_dt, dcB_dt};

This should reject the current parameter values immediately instead of waiting for the ODE to finish running, and it will give you debugging information on your parameters only when they’re bad. You should probably remove the print statements if you do this, as they would gum up the output stream.

I should also note that the syntax you’re using for the ODE functions (and for your arrays) are somewhat old; the array syntax at least will be deprecated somewhat soon, and the ODE syntax is much more cumbersome than the updated versions of the functions. I’d take a look at the ODE chapter in the Stan User’s Guide for information on the updated, more straightforward ODE interface. There are also examples of the new array syntax in the chapter.

Finally, it would be worth looking into vectorizing your likelihood, as likelihood functions are more efficient. If you define your observations and sigmas as matrix[3,T] objects, you could vectorize the likelihood as:

model {
    // define inits & such
    array[T] vector[3] y_sim = integrate_ode_rk45( ....... ); 
    matrix[3,T] y_sim_mat; // define a matrix version
    // Insert your priors here
        //...
    // Likelihood
    for(t in 1:T) {
        y_sim_mat[,t] = y_sim[t];
    }
    // vectorized likelihood, where y_obs and sigma are matrices defined in data block
    to_vector(y_obs) ~ normal(to_vector(y_sim_mat), to_vector(sigma));
}

Edit 2: My initial edit was inaccurate and has been removed.

1 Like
3

Hi @Christopher-Peterson,
Thank you very much for your response and suggestions for updating my code. If I am understanding correctly, the 'nan’s arise from the ODE dynamics pushing the variables into impossible space. I will attempt to implement the reject statement into the model.

I will have to dig in a bit more to the new syntax for ODE systems. I originally was using the updated syntax and had compiling issues (maybe from using rstan vs cmdstan to compile???).

Hopefully, I can update on my successes.

Thanks again.

4

That’s one possibility for the NaN’s. Adding the reject should let you know if that’s the issue, at least. Note: My initial edit about not being able to put the reject statement inside the ODE function was inaccurate, so that would probably be the easiest place to add it.

I’m not sure whether the latest version of rstan works with the new ODE syntax, but cmdstanr definitely does.

5

Hi @jimjunker1 , I am having a similar issue and was wondering if adding the reject statement into your model was successful for you?
I haven’t had any success, I have similar syntax as you but it appears to work fine with another model I have…
let me know if you might have any insight
Thanks

6

@morg_ana Apologies on the slow reply. Unfortunately, this did not immediately fix the issue and I have not had a chance to dive into it deeper. If I do, I will update with any progress here. Sorry it is not better news and best of luck.