Problems with bayesian ODEs parameter estimation

Questions
1

Hi all, I’m trying to use PyMC3 and specifically PyMC3.ode.DifferentialEquation to estimate some parameters of a ODE system.

I’m following the tutorial: https://docs.pymc.io/notebooks/ODE_API_introduction.html

My ODE system is called SEIR which is an epidemiological model, this system is normalized, which means that the sum of every state of the ODE is equal to 1, since every state is divided by the number of the population.

And it has 3 parameters I’m interested in: alpha, beta, and gamma.

I let alpha fixed = 0.2 and then I followed the next code to define my ODE system:

def SEIR(y, t, p):
"""
y[0] = S
y[1] = E
y[2] = I
y[3] = R


p[0] = Alpha
p[1] = Beta
p[2] = Gamma


St = S[-1] - (beta*S[-1]*I[-1])*dt
Et = E[-1] + (beta*S[-1]*I[-1] - alpha*E[-1])*dt
It = I[-1] + rho*(alpha*E[-1] - gamma*I[-1])*dt
Rt = gamma*I[-1]
"""



ds = -p[1]*y[0]*y[2]
de = p[1]*y[0]*y[2] - p[0]*y[1]
di = p[0]*y[1] - p[2]*y[2]
dr = p[2]*y[2]
return [ds, de, di, dr]

Then I simply followed the tutorial:

seir_model = DifferentialEquation(
func=SEIR,
times=times,
n_states=4,
n_theta=3,
)

And this is how I defined the probabilistic part of the model:

with pm.Model() as model:

sigma = pm.HalfCauchy('sigma',1, shape=4)   
p_gamma = pm.Bound(pm.Normal, lower=0, upper=1)("gamma",.5, 0.37 )  
R0 = pm.Bound(pm.Normal, lower=1, upper=5)('R0', 3.5, 1)
p_beta = pm.Deterministic('beta', gamma*R0)
alpha = 0.2


seir_curves = seir_model(y0=obs[0], theta=[alpha, p_beta, p_gamma])
Y = pm.Normal('Y', mu=seir_curves, sd=sigma, observed=obs)

prior = pm.sample_prior_predictive()
trace = pm.sample(2000, tune=500)
posterior_predictive = pm.sample_posterior_predictive(trace)
data = az.from_pymc3(trace=trace, prior = prior, posterior_predictive = posterior_predictive)

Regarding sigma, I have the following question since it is a normalized model, I should put as the parameter of the HalfCauchy something like:

(expected deviated cases) / Population

?

Also, the sampling process is super slow, so I don’t know if it has a relationship with how the model is specified.

Any suggestions are well-received!

Please, I need help with this.

2

I’m still having this speed issue. It’s super slow even in the toy example in the notebook.

3

Hi,
The performance issue is somewhat diffusely related to how we implemented the DifferentialEquation with theano. There are three options:

  • one can use a different sampler, that doesn’t use gradients (wastes performance on gradients though)
  • one can use a custom Op to wrap around scipy.integrate.odeint
  • @aseyboldt is working on sunode which is 10-50x faster than what we have now. The API is not yet stable, but there’s a LotkaVoltera ODE example in the README.

I strongly believe that option 3 is the way to go. Unfortunately, I don’t have enough time these days to be of much help.

1 Like