Structure of stan code block for deterministic ODE model

Modeling
,
1

Hi everyone

I am totally new to stan/rstan and I am considering fitting my compartmental, deterministic ODE model in stan (HMC).

  1. In which code block does the solver function (integrate_ode_rk45()) go? in the example in the user guide on this website, the integration happens in the model{} block, but every example I found online has the integration in transformed parameters{}.
  2. I have a quite complex model, which first runs an ODE model until the states reach equilibrium. This means I simulate for a prespecified time and then check the resulting states. The simulation is extended until some conditions are satisfied.
    These equilibrium states are then used as inits in a second ODE model. Is it possible to code this in stan? The fitting happens on one state variable in the second ODE model.

Thank you for your help!
Best, Fabienne

1 Like
2

Welcome!

  1. The main difference between model and transformed parameters is that transformed parameters are saved as part of the output. You do the integration in transformed parameters if you want to examine the resulting states outside of Stan. If you only care about the parameters that go in then the model block is fine.
  2. Does equilibrium mean the state converges to a single value? Or a limit cycle? If former, Stan’s algebraic solver may be more suitable than the ODE solver.
1 Like
3

Hi
Thanks for the fast reply. That’s very useful information. I assume running the solver in the transformed parameters block is slower due to the overhead of writing the data from the states in each step?

My model is a vector-borne epidemic model. The first ODE model is needed to initiate the vector population with its different life stages before disease introduction. The abundance at the different states depends on a carrying capacity, which is fitted, and on external temperature forcing. Once the population is in “equilibrium”, the abundance of all states is periodical (1 year). I cannot know what the equilibrium states are unless I simulate it and compare the values to the prior values a year ago. I don’t think this can be solved with an algebraic solver?

My idea was to implement the first model in the model block, save the output (states) and use them as inits for the second (disease) model. Since code is executed sequentially, I assume this works (in theory)?

TIA, Fabienne

4

Yes, it works in theory. And I agree it would be difficult (though probably not impossible) to find the equilibrium using the algebraic solver. It’s worth trying out the ODE model.