adapt {adapt}R Documentation

Adaptive Numerical Integration in 2–20 Dimensions

Description

Integrates a scalar function over a multidimensional rectangle, i.e., computes

integral[l .. u] functn(t) d^n(t)

where l =lower, u =upper and n =ndim. Infinite rectangles are not allowed, and ndim must be between 2 and 20.

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 R function 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

Details

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 R runs 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.

Value

A list of class "integration" with components

value

the estimated integral

relerr

the estimated relative error; < eps argument if the algorithm converged properly.

minpts

the actual number of function evaluations

ifail

an error indicator. If ifail is not equal to 0, the function warns the user of the error condition.

See Also

integrate

Examples

## 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("\np = ", p, "\n------\n")
  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))
}

[Package adapt version 1.0-4 Index]