# Simple ODE model drawing crazy values

OSX 10.11.6
Pystan 2.17.0.0 (I had compiler problems with Pystan 2.18)
Python 3.6.8, Anaconda
GCC 4.2.1 Compatible Clang 4.0.1

I’m trying to use Stan to fit a super-simple toy ODE model of one variable:

\frac{dx}{dt} = \alpha - \gamma_d*x

The Stan code looks like this, including a couple of diagnostic prints:

// Generative model for a very simple constitutive expression ODE -- x is
// produced at some constant rate and degraded/diluted proportionally to
// its concentration. Detection noise is lognormal around the true value.

// params[1] = alpha (production rate)
// params[2] = gamma_d (degradation/dilution rate)

functions {
real[] constitutive_production(real t,
real[] xs,
real[] params,
real[] x_r,
int[] x_i) {
real dx_dt[1];

dx_dt[1] = params[1] - xs[1]*params[2];
return dx_dt;
}
}

data {
real<lower = 0> init[1];
int<lower = 1> n_times;
int<lower = 1> n_params;
int<lower = 1> n_vars;
real<lower = 0> t0;
real<lower = 0> ts[n_times];
real<lower = 0> xs[n_times];
}

transformed data {
real x_r[0];
int x_i[0];
}

parameters {
real<lower=0> params[n_params];
/*
1: alpha
2: gamma_d
*/

}

transformed parameters{
real<lower = 0> x_hat[n_times, n_vars]; // Output from the ODE solver
real<lower = 0> x0[n_vars]; // Initial conditions

print("alpha: ", params[1]);
print("gamma_d: ", params[2]);

x0 = init;

x_hat = integrate_ode_rk45(constitutive_production, x0, t0, ts, params,
x_r, x_i);
}

model {
params[1] ~ normal(2, 2.5); // alpha
params[2] ~ normal(4, 2.5); // gamma_d

for (t in 1:n_times) {
xs[t] ~ normal(x_hat[t, 1], 0.1) T[0,];
}
}

generated quantities {
real<lower = 0> sim_x_hats[n_times, 1];

sim_x_hats = integrate_ode_rk45(constitutive_production, init, -0.1,
ts, params, x_r, x_i);
}


Data are generated with the following Python code (using Pystan, sorry Rstan users!):

init    = np.array([1e-4])
n_times = 50
t_min   = 1e-2
t_max   = 1
ts      = np.linspace(t_min, t_max, n_times)

real_alpha = 4
real_gamma = 2
real_sigma = 0.1
real_floor = 0.02

def ode_func(t, xs):
return np.array([real_alpha - xs[0]*real_gamma])

ode_sol = scipy.integrate.solve_ivp(ode_func, (t_min, t_max),
init, dense_output = True)

x_hats = ode_sol.sol(ts)[0]
# x_obs  = np.random.lognormal(np.log(x_hats + real_floor), real_sigma)
x_obs  = np.abs(np.random.normal(x_hats, real_sigma))


Running with:

simu_data ={"init": init,
"n_times": n_times,
"n_params": 2,
"n_times": len(ts),
"n_vars": len(init),
"t0": 0,
"ts": ts,
"xs": x_obs}

fit = model.sampling(data=simu_data, test_grad = False,
iter=100, warmup=0, chains=5,
seed=4832, )#, algorithm="Fixed_param")


When I run this, typical run-time output looks like the following, and I get no movement in any chains. It also takes a very long time to run, which I suppose isn’t surprising if it’s spending a ton of time rejecting proposals.

Start runtime message

Informational Message: The current Metropolis proposal is about to be rejected because of the following issue:
Exception: Max number of iterations exceeded (1000000). (in ‘unkown file name’ at line 54)

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

alpha: 7.3441
gamma_d: 4.54642

alpha: 7.3441
gamma_d: 4.54642

alpha: 7.83536e-117
gamma_d: 1.17744e+47

Informational Message: The current Metropolis proposal is about to be rejected because of the following issue:
Exception: Max number of iterations exceeded (1000000). (in ‘unkown file name’ at line 54)

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

alpha: 7.3441
gamma_d: 4.54642

alpha: 7.3441
gamma_d: 4.54642

alpha: 3.15723e-116
gamma_d: 8.56896e+46

Informational Message: The current Metropolis proposal is about to be rejected because of the following issue:
Exception: Max number of iterations exceeded (1000000). (in ‘unkown file name’ at line 54)

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

alpha: 7.3441
gamma_d: 4.54642

alpha: 7.3441
gamma_d: 4.54642
End runtime message

It looks like the sampler is veering off to crazy parameter values, priors be damned. It also looks like the two parameters are strongly anticorrelated, which doesn’t make sense – even if this model is underspecified, I would expect \alpha and \gamma_d to correlate, not anticorrelate. Any idea what’s tripping me up here?

Try using the ode function where you have control over the ode tolerances and the maxiter parameter. Also give the stiff ode solver using bdf integration a go .

Using integrate_ode_bdf instead of integrate_ode_rk45 causes a compile-time error:

FileNotFoundError Traceback (most recent call last)
/stan_utility.py in compile_model(filename, model_name, **kwargs)
204 try:
–> 205 sm = pickle.load(open(cache_fn, ‘rb’))
206 except:

FileNotFoundError: [Errno 2] No such file or directory: ‘cached-model-226ec904ac955d225c43d621eab9f88e.pkl’

During handling of the above exception, another exception occurred:

ImportError Traceback (most recent call last)
in ()
2
----> 3 model = stan_utility.compile_model(‘simplest_ode.stan’)#, extra_compile_args=extra_compile_args)

/Users/sclamons/Dropbox/Murray_Lab/software/stan/stan_utility.py in compile_model(filename, model_name, **kwargs)
206 except:
–> 207 sm = pystan.StanModel(model_code=model_code)
208 with open(cache_fn, ‘wb’) as f:
209 pickle.dump(sm, f)

/Users/sclamons/anaconda/lib/python3.6/site-packages/pystan/model.py in init(self, file, charset, model_name, model_code, stanc_ret, boost_lib, eigen_lib, verbose, obfuscate_model_name, extra_compile_args)
317 os.dup2(orig_stderr, sys.stderr.fileno())
318
–> 319 self.module = load_module(self.module_name, lib_dir)
320 self.module_filename = os.path.basename(self.module.file)
321 # once the module is in memory, we no longer need the file on disk

48 pyximport.install()
49 sys.path.append(module_path)
—> 50 return import(module_name)
51 else:
52 import imp

Referenced from: /var/folders/nr/rylr3yls5cddd1cybsbpqvjr0000gn/T/tmputagdd43/stanfit4anon_model_226ec904ac955d225c43d621eab9f88e_8158380586008065412.cpython-36m-darwin.so
Expected in: flat namespace
in /var/folders/nr/rylr3yls5cddd1cybsbpqvjr0000gn/T/tmputagdd43/stanfit4anon_model_226ec904ac955d225c43d621eab9f88e_8158380586008065412.cpython-36m-darwin.so

warmup=0

You might want to have some warmup done

Nailed it! I didn’t realize the burnin time is used for adaptation. Runs smooth as butter now, and correctly estimates my two parameters.

Thanks for the fix!

This equation is no-where near stiffness.

If the ODE is the production model instead of part of prototyping, one would simply use the analytical solution.