Adding 1d integrator to Stan


#1

I’ve created the issue and started to work on merging the 1D integrator into the language, but I think there may be things we need to fix first.

Here’s the issue and the integrator in stan-dev/math.

The doc is a little unclear on types and the type names don’t match in the tparam, but that’s not a big deal.

Issue 1

Shouldn’t there be a data argument for real parameters to the function that aren’t differentiated? And shouldn’t there be a data argument for integer parameters? As is, we’ll have to promote real constants and won’t be able to use integer constants from the program.

I think we want this signature for the function that gets integrated, where theta is a parameter vector, x_r and x_i are real and integer data, and y is the variable being integrated:

real
f(real y, vector theta, real[]  x_r, int[] x_i)

Then the function would look like this, where all the variables mean the same thing as above, and lb and ub anre the integration bounds and abs_tol and rel_tol are integration tolerances.

real
integrate_definite(F f, real lb, real ub, 
    vector theta, real[] x_r, int[] x_i,
    real abs_tol, real rel_tol)

Issue 2

Shouldn’t the vector type for theta be real[] instead? That’s what we did with the other integrators.

Issue 3

Was there a way to use this to do unbounded integrals with transforms or should I just document the straight-ahead definite integral usage for the manual?


#2

Issue 1 and 2

I tested it now and it works fine for real[] too. I can PR an additional unit test for this. I can also create and add a PR for the additional arguments you ask in Issue 1.

Issue 3

I think the integrator itself does the transformation to [0, 1], but as far as I tested, it didn’t work well with inf doubles as limits.

Note: there’s already a branch called feature/integrate-function-parser-doc at stan-dev/stan which already includes some attempts to add the integration to the parser and function generator.


#3

This integrator transforms to a common scale before estimating the integral, but those transforms don’t necessarily scale well in the limit that the bounds move out to infinity. Instead we can add new functions to handle (-inf, inf) and (0, inf) with proper transforms internally.


#4

That’d be great. It still leaves open the issue of which signatures to support in the Stan language.

Thanks—I can work from that branch, then. It shouldn’t be hard to add at all once we settle on signatures.


#5

One thing I keep meaning to point out about this is that we would only need the signature that @randommm has included if we allowed data binding through lambdas. The problem is that you can’t bind off parameters this way as those need to be identified for the integrator to deal with.


#6

How about having:

Data-only versions (easy to add the functions, don’t know about the parser though):

real
integrate_1d_tsc_simple(F f, real lb, real ub, 
    real abs_tol, real rel_tol)

real
integrate_1d_tsc_data_only(F f, real lb, real ub,
    vector_like x_r, int[] x_i,
    real abs_tol, real rel_tol)

A with derivative versions:

real
integrate_1d_tsc_derivative(F f, real lb, real ub, 
    vector_like param, vector_like x_r, int[] x_i,
    real abs_tol, real rel_tol)

real
integrate_1d_tsc_explicit_derivative(F f, F g, real lb, real ub, 
    vector_like param, vector_like x_r, int[] x_i,
    real abs_tol, real rel_tol)

vector_like can be either vector, real[] or a single real

Suffixes like data_only are not strictly necessary (or at least not that verbose), but might make it easier to identify each version.

About the transforms for unbounded integration, would the be usual logit/exp of Stan? I think it might worth addressing that specific problem after this one and the parser/doc to avoid delaying this too much (for now a explicit “no unbounded integration” in the manual).


#7

My preference would be to start with just a single function in Stan:

real integrate_definite(F f, real lb, real ub, 
    real[] param, real[] x_r, int[] x_i,
    real abs_tol, real rel_tol);

What does the function g compute (what’s it’s signature)? Is that for user-defined gradients to accelerate the integrator?


#8

Exactly.

OK!