Hi everyone, I want to use Stan to estimate the parameter in an ordinary differential equation. However, I don’t know how to input the time series of the boundary condition. It is not a “function”, but a time series (i.e., discrete data points).
For example, the ODE about a species’s concentration in a well-mixed lake is:
V*dC/dt = Load_{in} - Q_{out}*C - k*C
where V is the lake’s volume, C is the species’s concentration, Load_{in} is the species loading into the lake, Q_{out} is the outflow, k is the decay rate, t is time.
If the Load_{in} is a time-series data, e.g., daily load values, how can I write it into the Stan? (or any other software that can use this kind of boundary condition).
Hello @Tom8, if I understand you correctly, then the process Load_{in} in this model represents an input of your organism into the system. You observed the size of this input at some specific points in time, along with C, and you want to pass that as data, presumably to estimate the unknown parameters in the derivative function.
This seems analogous to the dose in pharmacokinetic models (which are also models of concentration of a well-mixed fluid). Typically in those applications the times of the dose are known, and the integration of the ODE system takes place between those times of dosing. At the end of each interdose period the ODE integration stops and the dosing is implemented as a change in the initial conditions for the ODE system for the next interdosing period. You just need to implement some indexing to match the sequence of doses, and times of dosing, with the desired output times from the ODE. I’d do this by embedding the ODE call in a for loop, along with the changes to the initial conditions for the next iteration. I don’t have an example to hand as I haven’t implemented a mutiple-dose system in Stan yet, though I suppose the principle is the same for any language. Perhaps there is some suitable code implemented in Torsten (https://discourse.mc-stan.org/t/bayesian-pharmacometrics-modeling-with-stan-and-torsten-preprint).
One point to consider is that you will need some substantive assumption about the time span over which Load_{in} is implemented. If it’s set up as ‘daily load’ then that sounds like a single point change (a ‘bolus dose’), which may not be realistic (depending on the relative time scale of events in this system).
I have more experience in structural dynamics/vibrations, so I wrote some code to illustrate my point. diff-eqs.html (878.5 KB)
It’s not a perfect solution because I have assumed you can write a function defining your forcing terms. One possible workaround that I may get to or not is modifying the user defined forcing function to accept an array of discretized forces and times, lookup the current time step, and return the value inside the ode_ call.