| 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 |
upper |
vector of at least length |
minpts |
the minimum number of function evaluations. |
maxpts |
the maximum number of function evaluations or
|
functn |
an R function which should take a single vector
argument and possibly some parameters and return the function value
at that point. |
eps |
the desired accuracy for the relative error. |
... |
other parameters to be passed to |
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; |
minpts |
the actual number of function evaluations |
ifail |
an error indicator. If |
## 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))
}