divonne {cubature}  R Documentation 
Divonne works by stratified sampling, where the partioning of the integration region is aided by methods from numerical optimization.
divonne(
f,
nComp = 1L,
lowerLimit,
upperLimit,
...,
relTol = 1e05,
absTol = 1e12,
minEval = 0L,
maxEval = 10^6,
flags = list(verbose = 0L, final = 1L, keep_state = 0L, level = 0L),
rngSeed = 0L,
nVec = 1L,
key1 = 47L,
key2 = 1L,
key3 = 1L,
maxPass = 5L,
border = 0,
maxChisq = 10,
minDeviation = 0.25,
xGiven = NULL,
nExtra = 0L,
peakFinder = NULL,
stateFile = NULL
)
f 
The function (integrand) to be integrated as in
This information might be useful if
the integrand takes long to compute and a sufficiently accurate
approximation of the integrand is available. The actual value
of the integral is only of minor importance in the partitioning
phase, which is instead much more dependent on the peak
structure of the integrand to find an appropriate
tessellation. An approximation which reproduces the peak
structure while leaving out the fine details might hence be a
perfectly viable and much faster substitute when

nComp 
The number of components of f, default 1, bears no relation to the dimension of the hypercube over which integration is performed. 
lowerLimit 
The lower limit of integration, a vector for hypercubes. 
upperLimit 
The upper limit of integration, a vector for hypercubes. 
... 
All other arguments passed to the function f. 
relTol 
The maximum tolerance, default 1e5. 
absTol 
the absolute tolerance, default 1e12. 
minEval 
the minimum number of function evaluations required 
maxEval 
The maximum number of function evaluations needed, default 10^6. Note that the actual number of function evaluations performed is only approximately guaranteed not to exceed this number. 
flags 
flags governing the integration. The list here is exhaustive to keep the documentation and invocation uniform, but not all flags may be used for a particular method as noted below. List components:

rngSeed 
seed, default 0, for the random number
generator. Note the articulation with 
nVec 
the number of vectorization points, default 1, but can be set to an integer > 1 for vectorization, for example, 1024 and the function f above needs to handle the vector of points appropriately. See vignette examples. 
key1 
integer that determines sampling in the partitioning
phase: 
key2 
integer that determines sampling in the final
integration phase: same as 
key3 
integer that sets the strategy for the refinement
phase: 
maxPass 
integer that controls the thoroughness of the
partitioning phase: The partitioning phase terminates when the
estimated total number of integrand evaluations (partitioning
plus final integration) does not decrease for 
border 
the relative width of the border of the integration
region. Points falling into the border region will not be
sampled directly, but will be extrapolated from two samples
from the interior. Use a nonzero 
maxChisq 
the maximum 
minDeviation 
a bound, given as the fraction of the requested
error of the entire integral, which determines whether it is
worthwhile further examining a region that failed the

xGiven 
a matrix ( 
nExtra 
the maximum number of extra points the peakfinder
subroutine will return. If 
peakFinder 
the peakfinder subroutine. This R function is
called whenever a region is up for subdivision and is supposed
to point out possible peaks lying in the region, thus acting as
the dynamic counterpart of the static list of points supplied
in 
stateFile 
the name of an external file. Vegas can store its entire internal state (i.e. all the information to resume an interrupted integration) in an external file. The state file is updated after every iteration. If, on a subsequent invocation, Vegas finds a file of the specified name, it loads the internal state and continues from the point it left off. Needless to say, using an existing state file with a different integrand generally leads to wrong results. Once the integration finishes successfully, i.e. the prescribed accuracy is attained, the state file is removed. This feature is useful mainly to define ‘checkpoints’ in longrunning integrations from which the calculation can be restarted. 
Divonne uses stratified sampling for variance reduction, that is, it partitions the integration region such that all subregions have an approximately equal value of a quantity called the spread (volume times halfrange).
See details in the documentation.
A list with components:
the actual number of subregions needed
the actual number of integrand evaluations needed
if zero, the desired accuracy was reached, if 1, dimension out of range, if 1, the accuracy goal was not met within the allowed maximum number of integrand evaluations.
vector of length nComp
; the integral of
integrand
over the hypercube
vector of
length nComp
; the presumed absolute error of
integral
vector of length nComp
;
the \chi^2
probability (not the
\chi^2
value itself!) that error
is not a
reliable estimate of the true integration error.
J. H. Friedman, M. H. Wright (1981) A nested partitioning procedure for numerical multiple integration. ACM Trans. Math. Software, 7(1), 7692.
J. H. Friedman, M. H. Wright (1981) User's guide for DIVONNE. SLAC Report CGTM193REV, CGTM193, Stanford University.
T. Hahn (2005) CUBAa library for multidimensional numerical integration. Computer Physics Communications, 168, 7895.
integrand < function(arg, phase) {
x < arg[1]
y < arg[2]
z < arg[3]
ff < sin(x)*cos(y)*exp(z);
return(ff)
}
divonne(integrand, relTol=1e3, absTol=1e12, lowerLimit = rep(0, 3), upperLimit = rep(1, 3),
flags=list(verbose = 2), key1= 47)
# Example with a peakfinder function
nDim < 3L
peakf < function(bounds, nMax) {
# print(bounds) # matrix (ndim,2)
x < matrix(0, ncol = nMax, nrow = nDim)
pas < 1 / (nMax  1)
# 1ier point
x[, 1] < rep(0, nDim)
# Les autres points
for (i in 2L:nMax) {
x[, i] < x[, (i  1)] + pas
}
x
} #end peakf
divonne(integrand, relTol=1e3, absTol=1e12,
lowerLimit = rep(0, 3), upperLimit = rep(1, 3),
flags=list(verbose = 2), peakFinder = peakf, nExtra = 4L)