Adaptive Gauss-Kronrod Quadrature for multiple limits

Integrates fun over the bounds [ lower, upper ] vectorized over lower and upper. Vectorized list structures of parameters can also be passed.

Usage

integrate_gk(
  fun,
  lower,
  upper,
  params = list(),
  .tolerance = .Machine$double.eps^0.25,
  .max_iter = 100L
)

Arguments

fun

A function to integrate. Must be vectorized and take one or two arguments, the first being points to evaluate at and the second (optionally) being parameters to apply. It must return a numeric vector the same length as its first input.

Currently, infinite bounds are not supported.

lower, upper

Integration bounds. Must have the same length.

params

Parameters to pass as a second argument to fun. The actual parameters must have the same length as the number of integrals to compute. Can be a possibly nested list structures containing numeric vectors. Alternatively, can be a matrix with the same number of rows as the number of integrals to compute.

.tolerance

Absolute element-wise tolerance.

.max_iter

Maximum number of iterations. The number of integration intervals will be at most length(lower) * .max_iter. Therefor the maximum number of function evaluations per integration interval will be 15 * .max_iter.

Value

A vector of integrals with the i-th entry containing an approximation of the integral of fun(t, pick_params_at(params, i)) dt over the interval lower[i] to upper[i]

Details

The integration error is estimated by the Gauss-Kronrod quadrature as the absolute difference between the 7-point quadrature and the 15-point quadrature. Integrals that did not converge will be bisected at the midpoint. The params object will be recursively subsetted on all numeric vectors with the same length as the number of observations.

Examples

# Argument recycling and parallel integration of two intervals
integrate_gk(sin, 0, c(pi, 2 * pi))
#> [1]  2.000000e+00 -3.141135e-16

dist <- dist_exponential()
integrate_gk(
  function(x, p) dist$density(x, with_params = p),
  lower = 0, upper = 1:10,
  params = list(rate = 1 / 1:10)
)
#>  [1] 0.6321206 0.6321206 0.6321206 0.6321206 0.6321206 0.6321206 0.6321206
#>  [8] 0.6321206 0.6321206 0.6321206
dist$probability(1:10, with_params = list(rate = 1 / 1:10))
#>  [1] 0.6321206 0.6321206 0.6321206 0.6321206 0.6321206 0.6321206 0.6321206
#>  [8] 0.6321206 0.6321206 0.6321206