adapt: Adaptive Numerical Integration in 2--20 Dimensions

Description

Integrates a scalar function over a multidimensional rectangle.

Usage

adapt(ndim, lower, upper, minpts = 100, maxpts = NULL, functn, 
    eps = 0.01, ...)

Arguments

ndim

the dimension of the integral, andi.e. number

lower

vector of at least length ndim of the lower bounds on the integral.

upper

vector of at least length ndim of the upper bounds on the integral.

minpts

the minimum number of function evaluations.

maxpts

the maximum number of function evaluations or NULL per default, see Details.

functn

an Rfunction which should take a single vector argument and possibly some parameters and return the function value at that point. functn must return a single numeric value.

eps

the desired accuracy for the relative error.

...

other parameters to be passed to functn

Value

A list of class "integration" with components
valuethe estimated integral
relerrthe estimated relative error; < eps argument if the algorithm converged properly.
minptsthe actual number of function evaluations
ifailan error indicator. If ifail is not equal to 0, the function warns the user of the error condition.

Details

The function computes $$\int_l^u \mbox{functn}(t)\,d^nt$$ where $l =$lower, $u =$upper and $n =$ndim. Infinite rectangles are not allowed, and ndim must be between 2 and 20. This is modified from Mike Meyer's S code. The functions just call A.C. Genz's fortran ADAPT subroutine to do all of the calculations. A work array is allocated within the C/Fortran code.

The Fortran function has been modified to use double precision, for compatibility with R. It only works in two or more dimensions; for one-dimensional integrals use the integrate function in the base package.

Setting maxpts to NULL asks the function to keep doubling maxpts (starting at max(minpts,500, r(ndim))) until the desired precision is achieved or Rruns out of memory. Note that the necessary number of evaluations typically grows exponentially with the dimension ndim, and the underlying code requires maxpts >= r(ndim) where $r(d) = 2^d + 2 d(d + 3) + 1$.

Examples

Run this code

## Example of  p - dimensional spherical normal distribution:
   ir2pi <- 1/sqrt(2*pi)
   fred <- function(z) { ir2pi^length(z) * exp(-0.5 * sum(z * z))}

   adapt(2, lo = c(-5, -5), up = c(5, 5), functn = fred)
   adapt(2, lo = c(-5, -5), up = c(5, 5), functn = fred, eps = 1e-4)
   adapt(2, lo = c(-5, -5), up = c(5, 5), functn = fred, eps = 1e-6)

## adapt "sees" function ~= constantly 0 --> wrong result
   adapt(2, lo = c(-9,-9), up = c(9,9), functn = fred)

## fix by using much finer initial grid:
   adapt(2, lo = c(-9,-9), up = c(9,9), functn = fred, min = 1000)
   adapt(2, lo = c(-9,-9), up = c(9,9), functn = fred, min = 1000, eps = 1e-6)
   i1 <- print(integrate(dnorm, -2, 2))$value

## True values for the following example:
   i1 ^ c(3, 5)
   for(p in c(3, 5)) {
     cat("p = ", p, "------
")
     f.lo <- rep(-2., p)
     f.up <- rep(+2., p)
     # not enough evaluations:
     print(adapt(p, lo=f.lo, up=f.up, max=100*p, functn = fred))
     # enough evaluations:
     print(adapt(p, lo=f.lo, up=f.up, max=10^p,  functn = fred))
     # no upper limit; p=3: 7465 points, ie 5 attempts (on an Athlon/gcc/g77):
     print(adapt(p, lo=f.lo, up=f.up, functn = fred, eps = 1e-5))
  }