Hi all,

I’m working on a toy example to estimate parameters for the Droop equations:

P = V*(P0 - P) - rho*(P / (Kp + P))*X;
Qp = rho*(P / (Kp + P)) - mu*(Qp - q0);

X = X

*mu*(1 - q0/Qp) - d*X;

I can solve the equations in STAN for a given set of parameters easily enough, and then add noise to those measurements (I’ve cross-checked the solver with scipy.integrate.odeint and it’s the same answer). The issue is then when I turn around and try to *estimate* the parameters, the model stalls and doesn’t fit well (and I often get “Exceeded maximum iteration” errors). Does anyone have any advice on what I’m doing wrong? Code below:

```
#%%
import numpy as np
import pystan as pyst
import matplotlib.pyplot as plt
from scipy.integrate import odeint
#%% Set up the parameters and initial conditions
V = 1
P0 = 2
rho = 0.7
Kp = 0.2
mu = 0.7
q0 = 0.3
d = 0.5
init = [5, 1, 2]
params = (V, P0, rho, Kp, mu, q0, d)
t = np.arange(0, 100, 1)
t0 = t[0]
ts = t[1:]
#%% STAN SOLVER
droop_stan_s = """
functions{
real[] droop(real t,
real[] y,
real[] theta,
real[] x_r,
int[] x_i){
real dydt[3];
dydt[1] = theta[1]*(theta[2] - y[1]) - theta[3]*(y[1] / (theta[4] + y[1]))*y[3];
dydt[2] = theta[3]*(y[1] / (theta[4] + y[1])) - theta[5]*(y[2] - theta[6]);
dydt[3] = y[3]*theta[5]*(1 - theta[6]/y[2]) - theta[7]*y[3];
return dydt;
}
}
data{
int<lower=1> T;
real t0;
real ts[T];
real y0[3];
real theta[7];
}
transformed data{
real x_r[0];
int x_i[0];
}
transformed parameters{
}
generated quantities{
real yhat[T,3];
yhat = integrate_ode_rk45(droop, y0, t0, ts, theta, x_r, x_i);
}
"""
droop_comp_s = pyst.StanModel(model_code=droop_stan_s)
#%% Simulate data
solver_dict = {'T': len(ts), 't0': t0, 'ts': ts, 'y0': init, 'theta': params}
sim_data = droop_comp_s.sampling(data=solver_dict,
algorithm='Fixed_param',
iter=1, chains=1)
yhat = sim_data.extract('yhat')['yhat']
plt.plot(ts, yhat[0,:,0], label='P')
plt.plot(ts, yhat[0,:,1], label='Qp')
plt.plot(ts, yhat[0,:,2], label='X')
plt.legend()
plt.show()
# add noise
yhat_obs = yhat + np.random.normal(0, 0.1, yhat.shape)
plt.plot(ts, yhat_obs[0,:,0], 'o', label = 'P')
plt.plot(ts, yhat_obs[0,:,1], 'o', label='Qp')
plt.plot(ts, yhat_obs[0,:,2], 'o', label='X')
plt.legend()
plt.show()
#%% STAN estimator
droop_stan_e = """
functions{
real[] droop(real t,
real[] y,
real[] theta,
real[] x_r,
int[] x_i){
real dydt[3];
dydt[1] = theta[1]*(theta[2] - y[1]) - theta[3]*(y[1] / (theta[4] + y[1]))*y[3];
dydt[2] = theta[3]*(y[1] / (theta[4] + y[1])) - theta[5]*(y[2] - theta[6]);
dydt[3] = y[3]*theta[5]*(1 - theta[6]/y[2]) - theta[7]*y[3];
return dydt;
}
}
data{
int<lower=1> T;
real y[T,3];
real t0;
real ts[T];
}
transformed data{
real x_r[0];
int x_i[0];
}
parameters{
real<lower=0> y0[3];
real theta[7];
vector<lower=0>[3] sigma_y;
}
transformed parameters{
}
model{
real yhat[T,3];
yhat = integrate_ode_rk45(droop, y0, t0, ts, theta, x_r, x_i);
for (t in 1:T){
y[t] ~ normal(yhat[t], sigma_y);
}
sigma_y ~ cauchy(0, 2.5);
theta ~ normal(0,1);
y0 ~ normal(0,1);
}
"""
droop_comp_e = pyst.StanModel(model_code=droop_stan_e)
#%%
estimate_dict = {'T': len(ts), 't0': t0, 'ts': ts, 'y': yhat_obs[0]}
droop_output = droop_comp_e.sampling(data=estimate_dict,
iter=10)
Informational Message: The current Metropolis proposal is about to be rejected because of the following issue:
Exception: Max number of iterations exceeded (1000000). (in 'unknown file name' at line 34)
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.
Informational Message: The current Metropolis proposal is about to be rejected because of the following issue:
Exception: Max number of iterations exceeded (1000000). (in 'unknown file name' at line 34)
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.
```

Also, can someone tell me what x_r and x_i are? I’ve read the docs, but I’m confused as to what they are what they are set to a 0-array (and sometimes I set x_i set to a 1-array).

Thanks, I appreciate the help!