Numerical Integrals

Hi Team,

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.

In the ONR proposal, we look at QUADPACK (see and In particular, the routines of interest are:

  • QAGI for integration over infinite intervals
  • QAGP for integration of functions with singularities and discontinuities

These routines use adaptive Gauss-Kronrod quadrature methods. The source code is in Fortran (, with a public domain license. The GNU Scientific Library reimplemented the QUADPACK routines in C ( under the GNU (general public) license.

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:

  1. You can always integrate over a large range. The price to pay for this should not be too costly as the integrators work adaptivley.

  2. 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.


integrate_ode.stan (481 Bytes)
integration.R (1.8 KB)

@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?

Glad i was able to help… what do you mean with not being able to call integrate_ode_bdf? In a stan model?

Maybe i should really collect those ode tricks in a document…

Yes, basically the model does not compile on Pystan, it just gives:

ImportError: /tmp/tmpfps6cgju/ undefined symbol: CVodeSetUserData

It works fine for integrate_ode_rk45, though.

Yup. PyStan doesn’t yet ship the library.

I have just caught up with this thread and the impressive code. One thing I don’t understand is the purpose of the 1e-5 in the start time.

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).

Of course. I see that now. It is a very fast integrator but not user-friendly.

…better than being stuck, I guess… but as I said, rescue is on its way.