Starting a discussion about doing numerical integration in Stan. For continuous functions with finite bounds, we can use ODE integrators, but we’ll need another tool for infinite intervals and functions with singularities and discontinuities.
Can we have an example where the ODE integrator fails?
I am not opposing this, but I guess that we get very far wit the ODE integrators:
You can always integrate over a large range. The price to pay for this should not be too costly as the integrators work adaptivley.
As long as you have a clue where the singularities are you can simply request a value to be output at the singularity. This will make the ODE integrator stop at this point and do a bit extra work needed.
That said, the ODE integrator is clearly a compromise and will most likely be very slow in comparison (there is a reason for QUADPACK being there) - but it works now already. Another approach is to simply sample that integral. Mike described that approach once.
Inspired by the tanh-sinh stuff pointed out from Ben, I figured a way to abuse our ODE integrator to do indefinite integration. Essentially by transforming from -Inf - +Inf to 0 - 1 with the inv_logit function. Probably there are smarter transforms.
I also added a R version of the sinh-tanh integration scheme. This R version can be directly translated to Stan (vanilla Stan functions code). Performance won’t be great, but at least it runs…no idea if the integrator is faster or this scheme. The speed will depend on how clever one can take advantage of vectorization to avoid a large AD stack. The only drawback is no adaptive error control, so one would have to sample from the prior and fix m and beforehand, but that is doable I think.
@wds15: Thanks for sharing the code. This seems to work pretty well!
Tbh, I’m not sure what kind of speed up we would get with numerical integration and @billg is probably the best person to answer your questions about the kind of problems we may run into if we stick the ODE integrator.
Great! I didn’t knew it was possible to use the ODE solver for numerical integration (yeah, in the end I know very little about the subject). After some attempts I was able to make it work, for now it’s way better than using my own code since I don’t need to ask users for a patched version of pystan.
For some odd reason, though, it seems I’m unable to call integrate_ode_bdf, receiving some kind undefined symbol. @ariddell is this a known issue?
Oh… glad you like it. I think the 1E-5 keeps it away from the boundary 0 and 1 where things get singular… but I don’t have the files open to confirm. If you paste the original code here, I can look again.
Using the ODE integrator as an integration function is tricky and often overlooked, yes. We will have an integrate_1d function (or similarly called) very soon in Stan (if it is not already in 2.18; I haven’t checked).