Hierarchical Lotka Volterra

With more than two populations, especially for long integration times, the ODE integration is unlikely to properly calculate the chaotic trajectories and it will be hard to say much of anything that depends on the precise details of the solution at those later times. You can try increasing the absolute and relative precision of the integration, but you’re fighting exponential errors so it’s a losing game. Chaotic systems are hard.

Thanks for the feedback - much appreciated - I will try a Lorenz system to see what happens - maybe this sort of problem is better suited to a Sequential Monte Carlo approach?

I’ve just pushed an update of ctsem to cran that allows this sort of system specified in stan as a nonlinear stochastic differential equation. HMC is very slow but optimizing and then importance sampling (routine included) seems to work well. So far no documentation and limited testing for the nonlinear aspect (integration using unscented kalman filters / runge kutta 4), but I’d be happy to chat / pass on examples if anyone wants to try it out…

2 Likes

It’s the ODE, not the MCMC that’s causing the problem here. So SMC isn’t going to help if you still have to solve a tricky ODE.

What I meant was treating the parameter as state. I’ve written an example of estimating a parameter for a chaotic Lorenz system using libbi:

http://idontgetoutmuch.org/estimating-chaotic.htm

Note that parameter is updated as log Brownian motion.

sub transition(delta = h) {
  w ~ normal (0.0, sqrt(h))
  ode(h = h, atoler = delta_abs, rtoler = delta_rel, alg =  'RK4(3)') {
    dX/dt = exp(ln_alpha) * (Y - X)
    dY/dt = X * (rho - Z) - Y
    dZ/dt = X * Y - beta * Z
    dln_alpha/dt = -sigma * sigma / 2 - sigma * w / h
  }
}

I am not sure how to do this in Stan but will go and study the manual.

EDIT: The Jupyter notebook is here

Is this available to try out?