Fitting an ODE time course in Stan

[Hua Xie]

I recently find out about several changes to a project I am working on in graduate school, and rather than fitting a vector of stable ODE solutions, I am now trying to fit time courses for which some input data is changing at every time step in my equations, i.e. for a dx/dt = I + g(x), I is changing at every dt. Actual model is a lot more complex than that, of course.

I am figuring out the best way of doing this, as I will not be able to use integrate_ode for my purposes because of the transient input data. I am debating between coding my entire script in Python/Julia with a simple Metropolis-Hastings MCMC algorithm for parameter fitting, trying out PyMC, or doing a simple Euler in Stan and still taking advantage of Stan's HMC. I just wanted to see if people had some success stories/example code of coding their own simple ODE solvers in Stan.

Thanks

[Michael Betancourt]

Having time-dependent forcing terms, dx/dt = I(t) + g(x) is in general

not an issue for integrate_ode.  If you know a functional form for

I(t) then you just add it to the ODE functor specification.  If I(t) is

piecewise constant then you can recreate that using if/then statements

in the ODE functor specification.

Alternatively you can call integrate_ode only over the periods for

which I(t) is constant, instead of trying to integrate over many points

at once.

You would have to implement the ODE this way no matter where

you build your analysis, so moving to Python or Julia (or, gulp,

trying to get away with an Euler integrator) isn’t going to make

that any easier than building it in Stan.  Plus you’d immediately

loose the scalable MCMC and be introducing even worse problems.

--
You received this message because you are subscribed to the Google Groups "Stan users mailing list" group.
To unsubscribe from this group and stop receiving emails from it, send an email to stan-users+...@googlegroups.com.
To post to this group, send email to stan-...@googlegroups.com.
For more options, visit https://groups.google.com/d/optout.

[Hua Xie]

Hi Michael,

Thanks for your help. I have started out getting the most basic of models working in Stan to learn it, but admittedly am still confused and have not been able to find much in the way of examples with forcing. I actually have multiple forcing terms in my actual mode, some piecewise linear, others for which I have a functional form. Here is a simplified template of the data and model I am working with:

int N //Where N = 40. My data is sampled on 50 days over a 500 day interval
real C[N] //Time course data I am trying to fit
real I[N] //Data for input forcing piecewise constant function that changes with every new N.
real Z[T] //Forcing for which I have functional form. Let's say Z[T] = sin(T).

dX/dT = I - theta[1] * X + theta[2] * Y / Z * (1 + X + Y)
dY/dT = -theta[3] * Y + theta[4] * X / Z * (1 + X + Y)

Chat[N] = (1 - theta[1]) * X[N]

C[N] ~ normal(Chat[N],sigma)

If you can offer me some guidance and clarification in terms of Stan pseudocode, I would really appreciate it. What would my ODE function block and integrate_ode statement look like with a piecewise constant function? Is this where I would use the x_r variable?

Thank you.

[Michael Betancourt]

real[] onecomp_mm_infusion(real t, real[] y, real[] theta, real[] x_r, int[] x_i) {
real dydt[1];

if (t < x_r[2])
dydt[1] <- - theta[2] * y[1] / ( theta[1] * (theta[3] + y[1]) ) + x_r[1] / (theta[1] * x_r[2]);
else
dydt[1] <- - theta[2] * y[1] / ( theta[1] * (theta[3] + y[1]) );

return dydt;

[Sebastian Weber]

Just as an additional comment to such ODEs: The ODE RHS function is discontinuous at t=x_r[2] as there is a sudden change. The integrated function itself is still continuous, but the ODE RHS is not. The latter is somewhat problematic for any integrator (ODE integrator always assume continuity of the integrator function and the integrated RHS) and one should work around it by stopping and restarting the integration there. So we want to integrate from t0 < tmid = x_r[2] < tend such a function with initial y0 at t0, then what is ideal to do is to run two commands of integrate_ode as:

part1 = integrate_ode(ode_rhs, t0, y0, ts1, ...); with ts1 = (... tmid), tmid is the n1'th time-point

then we do

part2 = integrate_ode(ode_rhs, tmid, part1[n1], ts2, ...); with ts2 = (one-after-tmid, ... tend)

That's it. Now the ODE RHS is always continuous when given to the integrator. Not doing this will likely work, but the solution won't be precise around the discontinuity and the integrator will be inefficient in that region as the step-size will be decreased to very small values.

Sebastian

[Hua Xie]

Thanks for the heads up on this. So I am not sure I will be able to do that particular method of starting and stopping the integrator with individual parts outside of the integrate_ode function, because I will have an x_r array with tens of thousands of changes (I feel that there are some problems with both the models and the data, but it's just what I have to work with, for now). Definitely anticipating inefficiencies, if having a step-wise forcing function with tens of thousands of discontinuities will even work at all.

[Hua Xie]

One more question --
So the ODE equations are on an hourly scale. One of the input terms changes every hour. My original plan was to have something like
dydt <- x_r[iceil(t)] - etc. However, Stan does not seem to have a ceiling/floor function that returns integers, last I checked. Is there an alternative to using a ceiling/floor function that comes to mind?

Thanks

[Sebastian Weber]

No, there is no ceiling function in Stan. Given your forcing function is that flexible and changing so often... maybe replace that with a spline or something else which is a continuous approximation to that forcing function?

The problem you describe sounds like you will have to think about how to make it easier to fit, i.e. approximate somehow and simplify.

Sebastian

[Hua Xie]

Yes, I would personally love to replace it with a continuous approximation of the forcing function. However, my supervisor still wants to fit the discontinuous data, and have been unable to convince him otherwise. I will continue discussing that with him -- in the mean time, is there any way to code what I described?

Thanks again for your patient help, Sebastian.

--
You received this message because you are subscribed to a topic in the Google Groups "Stan users mailing list" group.
To unsubscribe from this topic, visit https://groups.google.com/d/topic/stan-users/huPXv3NWKEU/unsubscribe.
To unsubscribe from this group and all its topics, send an email to stan-users+...@googlegroups.com.