11.4 1D integrator
Stan provides a built-in mechanism to perform 1D integration of a function via quadrature methods.
It operates similarly to the algebraic solver and the ordinary differential equations solver in that it allows as an argument a function.
Like both of those utilities, some of the arguments are limited to data only expressions. These expressions must not contain variables other than those declared in the data or transformed data blocks.
11.4.1 Specifying an integrand as a function
Performing a 1D integration requires the integrand to be specified somehow. This is done by defining a function in the Stan functions block with the special signature:
real integrand(real x, real xc, array real theta, array real x_r, array int x_i)
The function should return the value of the integrand evaluated at the point x.
The argument of this function are:
x, the independent variable being integrated over
xc, a high precision version of the distance from x to the nearest endpoint in a definite integral (for more into see section Precision Loss).
theta, parameter values used to evaluate the integral
x_r, data values used to evaluate the integral
x_i, integer data used to evaluate the integral
Like algebraic solver and the differential equations solver, the 1D
integrator separates parameter values,
theta, from data values,
11.4.2 Call to the 1D integrator
(function integrand, real a, real b, array real theta, array real x_r, array int x_i)
Integrates the integrand from a to b.
Available since 2.23
(function integrand, real a, real b, array real theta, array real x_r, array int x_i, real relative_tolerance)
Integrates the integrand from a to b with the given relative tolerance.
Available since 2.23
184.108.40.206 Arguments to the 1D integrator
The arguments to the 1D integrator are as follows:
integrand: function literal referring to a function specifying the integrand with signature
(real, real, array real, array real, array int):realThe arguments represent
- where integrand is evaluated,
- distance from evaluation point to integration limit for definite integrals,
- real data
- integer data, and the return value is the integrand evaluated at the given point,
a: left limit of integration, may be negative infinity, type
b: right limit of integration, may be positive infinity, type
theta: parameters only, type
x_r: real data only, type
x_i: integer data only, type
relative_tolerance argument can optionally be provided for more control over the algorithm:
relative_tolerance: relative tolerance for the 1d integrator, type
real, data only.
220.127.116.11 Return value
The return value for the 1D integrator is a
real, the value of the integral.
18.104.22.168 Zero-crossing integrals
For numeric stability, integrals on the (possibly infinite) interval \((a, b)\) that cross zero are split into two integrals, one from \((a, 0)\) and one from \((0, b)\). Each integral is separately integrated to the given
22.214.171.124 Precision loss near limits of integration in definite integrals
When integrating certain definite integrals, there can be significant precision loss in evaluating the integrand near the endpoints. This has to do with the breakdown in precision of double precision floating point values when adding or subtracting a small number from a number much larger than it in magnitude (for instance,
1.0 - x).
xc (as passed to the integrand) is a high-precision version of the distance between
x and the definite integral endpoints and can be used to address this issue. More information (and an example where this is useful) is given in the User’s Guide. For zero crossing integrals,
xc will be a high precision version of the distance to the endpoints of the two smaller integrals. For any integral with an endpoint at negative infinity or positive infinity,
xc is set to
126.96.36.199 Algorithmic details
Internally the 1D integrator uses the double-exponential methods in the Boost 1D quadrature library. Boost in turn makes use of quadrature methods developed in (Takahasi and Mori 1974), (Mori 1978), (Bailey, Jeyabalan, and Li 2005), and (Tanaka et al. 2009).
The gradients of the integral are computed in accordance with the Leibniz integral rule. Gradients of the integrand are computed internally with Stan’s automatic differentiation.