don’t have a lot of time… the fundamental issue is that the stan ode solver only gurantees you positivity up to an accuracy of the absolute tolerance at which you run the ODE integrator (the default is 1E-9 or so). Thus, the integrator will return negative values in that order. You have to solutions here:
-
easy: just add 1E-6 or so to the ode solution (so log(y_ode + 1E-6))… the 1E-6 needs to be reasonably above the absolute error at which you run the integrator… this is an easy solution, but more of a hack
-
better, but also has issues: right now you integrate dy/dt = f(t,y) and you get y(t) which is positive in your case… you can just make the integrator return log(y(t)) directly by integrating 1/y * dy/dt = d(log(y))/dt = f(y,t)/exp(log(y(t)))
the second approach is more elegant and more stable if you y does only vary smoothly in log space