| adapt {adapt} | R Documentation |
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.
adapt(ndim, lower, upper, minpts = 100, maxpts = NULL, functn, eps = 0.01, ...)
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 |
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.
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. |
## 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))
}