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?