ODE Model Not Converging

Modeling
1

Hi all,

I am trying to create a Stan model to perform parameter estimation of a basic SIR model using simulated data. However, when I run the model, I keep getting the following Informational Message:

‘The current Metropolis proposal is about to be rejected because of the following issue:
Exception: Max number of iterations exceeded (1000)’

The model is runs far more slowly than I feel it should, but does complete. However, it’s estimates are way off, and it finishes with a warning that 82.5% of the iterations ‘ended with a divergence’. I have tried increasing adapt_delta to 0.9, and reparameterizing my model but this hasn’t helped. Does anyone have any advice?

Here is my model:

import pystan
import numpy as np
import pandas as pd

model = """
functions{
real[] dz_dt(real t, real[] z, real[] theta, 
              real[] x_r, int[] x_i){
    real N = 1000;
    real s = z[1];
    real i = z[2];
    real r = z[3];
        
    real alpha = theta[1];
    real beta = theta[2];
        
    real ds_dt = -alpha * s * i * N;
    real di_dt = alpha * s * i * N- beta * i;
    real dr_dt = beta * i;
    return {ds_dt, di_dt, dr_dt};
}
}

              
data{
     int<lower = 0> M;          //Number of measurements
     real ts[M];                //Measurement times > 0
     real<lower = 0> y_init[3];  //Initial measured population
     real<lower = 0> z_init[3];  //True initial population
     real<lower = 0> y[M, 3];    //Measured population at measurement times
}

parameters{
    real<lower = 0> theta[2];   //theta = {alpha, beta}
    real<lower = 0> sigma;   //error scale
}

transformed parameters{
    real z[M, 3]
    = integrate_ode_rk45(dz_dt, z_init, 0.0, ts, theta,
                         rep_array(0.0,0), rep_array(0,0),
                         1e-6, 1e-5, 1e3);
}

model{
     theta[{1,2}] ~ gamma(1.5, 0.1);
     sigma ~ gamma(2.3, 0.01);
     for (k in 1:3) {
             y_init[k] ~ normal(z_init[k], sigma);
             y[ ,k] ~ normal(z[, k], sigma);
     }
}
"""

sm = pystan.StanModel(model_code = model)

num_meas = 100 #Number of measurements
M = num_meas - 1 #Number of measurements minus initial condition
t = np.arange(1,M+1)

alpha_true = 0.001
beta_true = 0.09

initial_inf = 1
N = 1000

y0 = [(N - initial_inf)/N ,initial_inf/N,0]

s,i,r= np.zeros(num_meas), np.zeros(num_meas), np.zeros(num_meas)
s_noise = np.zeros(num_meas)
i_noise = np.zeros(num_meas)
r_noise = np.zeros(num_meas)
tot = np.zeros(num_meas)
s[0], i[0], r[0] = y0

s_check, i_check, r_check = np.zeros(num_meas),np.zeros(num_meas),np.zeros(num_meas)
# s[1] = s[0] - alpha_true * s[0] * i[0] * N
# i[1] = i[0] + alpha_true * s[0] * i[0] * N - beta_true * i[0]
# r[1] = r[0] + beta_true * i[0]

for x in range(M):
    s[x+1] = s[x] - alpha_true * s[x] * i[x] * N
    i[x+1] = i[x] + alpha_true * s[x] * i[x] * N - beta_true * i[x]
    r[x+1] = r[x] + beta_true * i[x]
    
for x in range(M+1):
    s_noise[x] = min(max(0.0,  s[x]+np.random.normal(0,0.015)),1) 
    i_noise[x] = max(0.0, i[x]+np.random.normal(0,0.015))
    if s_noise[x] + i_noise[x] > 1:
        i_noise[x] = 1 - s_noise[x]
    r_noise[x] = 1 - s_noise[x] - i_noise[x]
    tot[x] = s_noise[x] + i_noise[x] + r_noise[x]

y0 = [(N - initial_inf)/N ,initial_inf/N,0]

data = np.array([s_noise[1:num_meas],i_noise[1:num_meas],r_noise[1:num_meas]])
data_tr = data.transpose()

stan_data = {'M': M, 'ts':t, 'y_init': y0,'z_init': y0, 'y':data_tr}

fit = sm.sampling(data = stan_data, chains = 2, iter = 1000, n_jobs=1, control = dict(adapt_delta=0.9))

Here is a copy-paste of my traceback:

If this warning occurs sporadically, such as for highly constrained variable typ
es like covariance matrices, then the sampler is fine,
but if this warning occurs often then your model may be either severely ill-cond
itioned or misspecified.

Informational Message: The current Metropolis proposal is about to be rejected b
ecause of the following issue:
Exception: Max number of iterations exceeded (1000). (in ‘unknown file name’ at
line 35)

If this warning occurs sporadically, such as for highly constrained variable typ
es like covariance matrices, then the sampler is fine,
but if this warning occurs often then your model may be either severely ill-cond
itioned or misspecified.

Informational Message: The current Metropolis proposal is about to be rejected b
ecause of the following issue:
Exception: Max number of iterations exceeded (1000). (in ‘unknown file name’ at
line 35)

If this warning occurs sporadically, such as for highly constrained variable typ
es like covariance matrices, then the sampler is fine,
but if this warning occurs often then your model may be either severely ill-cond
itioned or misspecified.

Informational Message: The current Metropolis proposal is about to be rejected b
ecause of the following issue:
Exception: Max number of iterations exceeded (1000). (in ‘unknown file name’ at
line 35)

If this warning occurs sporadically, such as for highly constrained variable typ
es like covariance matrices, then the sampler is fine,
but if this warning occurs often then your model may be either severely ill-cond
itioned or misspecified.

Informational Message: The current Metropolis proposal is about to be rejected b
ecause of the following issue:
Exception: Max number of iterations exceeded (1000). (in ‘unknown file name’ at
line 35)

If this warning occurs sporadically, such as for highly constrained variable typ
es like covariance matrices, then the sampler is fine,
but if this warning occurs often then your model may be either severely ill-cond
itioned or misspecified.

Informational Message: The current Metropolis proposal is about to be rejected b
ecause of the following issue:
Exception: Max number of iterations exceeded (1000). (in ‘unknown file name’ at
line 35)

If this warning occurs sporadically, such as for highly constrained variable typ
es like covariance matrices, then the sampler is fine,
but if this warning occurs often then your model may be either severely ill-cond
itioned or misspecified.

Informational Message: The current Metropolis proposal is about to be rejected b
ecause of the following issue:
Exception: Max number of iterations exceeded (1000). (in ‘unknown file name’ at
line 35)

If this warning occurs sporadically, such as for highly constrained variable typ
es like covariance matrices, then the sampler is fine,
but if this warning occurs often then your model may be either severely ill-cond
itioned or misspecified.

Informational Message: The current Metropolis proposal is about to be rejected b
ecause of the following issue:
Exception: Max number of iterations exceeded (1000). (in ‘unknown file name’ at
line 35)

If this warning occurs sporadically, such as for highly constrained variable typ
es like covariance matrices, then the sampler is fine,
but if this warning occurs often then your model may be either severely ill-cond
itioned or misspecified.

Informational Message: The current Metropolis proposal is about to be rejected b
ecause of the following issue:
Exception: Max number of iterations exceeded (1000). (in ‘unknown file name’ at
line 35)

If this warning occurs sporadically, such as for highly constrained variable typ
es like covariance matrices, then the sampler is fine,
but if this warning occurs often then your model may be either severely ill-cond
itioned or misspecified.

Informational Message: The current Metropolis proposal is about to be rejected b
ecause of the following issue:
Exception: Max number of iterations exceeded (1000). (in ‘unknown file name’ at
line 35)

If this warning occurs sporadically, such as for highly constrained variable typ
es like covariance matrices, then the sampler is fine,
but if this warning occurs often then your model may be either severely ill-cond
itioned or misspecified.

Informational Message: The current Metropolis proposal is about to be rejected b
ecause of the following issue:
Exception: Max number of iterations exceeded (1000). (in ‘unknown file name’ at
line 35)

If this warning occurs sporadically, such as for highly constrained variable typ
es like covariance matrices, then the sampler is fine,
but if this warning occurs often then your model may be either severely ill-cond
itioned or misspecified.

Informational Message: The current Metropolis proposal is about to be rejected b
ecause of the following issue:
Exception: Max number of iterations exceeded (1000). (in ‘unknown file name’ at
line 35)

If this warning occurs sporadically, such as for highly constrained variable typ
es like covariance matrices, then the sampler is fine,
but if this warning occurs often then your model may be either severely ill-cond
itioned or misspecified.

Informational Message: The current Metropolis proposal is about to be rejected b
ecause of the following issue:
Exception: Max number of iterations exceeded (1000). (in ‘unknown file name’ at
line 35)

If this warning occurs sporadically, such as for highly constrained variable typ
es like covariance matrices, then the sampler is fine,
but if this warning occurs often then your model may be either severely ill-cond
itioned or misspecified.

Informational Message: The current Metropolis proposal is about to be rejected b
ecause of the following issue:
Exception: Max number of iterations exceeded (1000). (in ‘unknown file name’ at
line 35)

If this warning occurs sporadically, such as for highly constrained variable typ
es like covariance matrices, then the sampler is fine,
but if this warning occurs often then your model may be either severely ill-cond
itioned or misspecified.

Informational Message: The current Metropolis proposal is about to be rejected b
ecause of the following issue:
Exception: Max number of iterations exceeded (1000). (in ‘unknown file name’ at
line 35)

If this warning occurs sporadically, such as for highly constrained variable typ
es like covariance matrices, then the sampler is fine,
but if this warning occurs often then your model may be either severely ill-cond
itioned or misspecified.

Informational Message: The current Metropolis proposal is about to be rejected b
ecause of the following issue:
Exception: Max number of iterations exceeded (1000). (in ‘unknown file name’ at
line 35)

If this warning occurs sporadically, such as for highly constrained variable typ
es like covariance matrices, then the sampler is fine,
but if this warning occurs often then your model may be either severely ill-cond
itioned or misspecified.

Informational Message: The current Metropolis proposal is about to be rejected b
ecause of the following issue:
Exception: Max number of iterations exceeded (1000). (in ‘unknown file name’ at
line 35)

If this warning occurs sporadically, such as for highly constrained variable typ
es like covariance matrices, then the sampler is fine,
but if this warning occurs often then your model may be either severely ill-cond
itioned or misspecified.

Iteration: 1000 / 1000 [100%] (Sampling)

Elapsed Time: 730.557 seconds (Warm-up)
499.38 seconds (Sampling)
1229.94 seconds (Total)

WARNING:pystan:825 of 1000 iterations ended with a divergence (82.5 %).
WARNING:pystan:Try running with adapt_delta larger than 0.9 to remove the diverg
ences.
Exectuion time 2773.0 secs
WARNING:pystan:Deprecation warning. PyStan plotting deprecated, use ArviZ librar
y (Python 3.5+). pip install arviz; arviz.plot_trace(fit))

2

Can you make plots of the ODE solution at various points in parameter space and see if they resemble what you expect?

1000 doesn’t seem like a large maximum number of timesteps. The default is 1e6.

You could also try the implicit solver (integrate_ode_bdf).

3

You have a scale problem: \alpha & \beta in your ODE system are in the same scale by both being sampled from gamma distribution, but N\gg \beta blows up the rate for s and i term. In other words, s & i change much faster than r, which causes the system to be stiff. As @bbbales2 said, try bdf integrator. Like he said, I wouldn’t imagine SIR model being terribly stiff except at t=0 so rk45 should also work when max time step is increased.

4

Thanks both for the advice. I redefined the model only existed in terms of s & i (I’d hoped that would fix the scale problem) but sadly did not help. I also tried switching to bdf integrator, but now my STAN model won’t compile. Advice very much appreciated nonetheless.

5

PyStan has trouble with the bdf integrator since this is non-header only. There are some instructions flying around to make this work. The easier route is to use CmdStanPy, I think - or another interface which supports the bdf integrator (CmdStan, RStan, CmdStanR, CmdStanPy).

Still, the scale mismatch topic is really important. Did you make plots of the system; maybe for the prior medians and some draws from the prior? That should help guide you in setting this up right.