# divonne: Integration by Stratified Sampling for Variance Reduction In R2Cuba: Multidimensional Numerical Integration

## Description

Divonne works by stratified sampling, where the partioning of the integration region is aided by methods from numerical optimization.

## Usage

 ```1 2 3 4 5 6 7 8 9``` ```divonne(ndim, ncomp, integrand, ..., lower=rep(0,ndim), upper=rep(1,ndim), rel.tol= 0.001, abs.tol = 0, flags=list(verbose=1, final=1, pseudo.random=0, mersenne.seed=NULL), min.eval=0, max.eval=50000, key1=47, key2=1, key3=1, max.pass=5, border=0, max.chisq=10, min.deviation=0.25, xgiven=NULL, nextra=0, peakfinder=NULL) ```

## Arguments

 `ndim` same as `cuhre` ` ncomp` same as `cuhre` ` integrand` same as `cuhre`. But, here, the input argument `phw` indicates the integration phase: `0`: sampling of the points in `xgiven`, `1`: partitioning phase, `2`: final integration phase, `3`: refinement phase. 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 integrand 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 `phw < 2`. In all other instances, `phw` can be ignored. ` ...` same as `cuhre` `lower` same as `cuhre` `upper` same as `cuhre` `rel.tol` same as `cuhre` `abs.tol` same as `cuhre` `flags` same as `cuhre` `pseudo.random` and `mersenne.seed` are only taken into account when the argument `key1` is negative. `min.eval` same as `cuhre` `max.eval` same as `cuhre` `key1` integer that determines sampling in the partitioning phase: `key1 = 7, 9, 11, 13` selects the cubature rule of degree `key1`. Note that the degree-11 rule is available only in 3 dimensions, the degree-13 rule only in 2 dimensions. For other values of `key1`, a quasi-random sample of \code{n=|key1|} points is used, where the sign of `key1` determines the type of sample, `key1 = 0`, use the default rule. `key1 > 0`, use a Korobov quasi-random sample, `key1 < 0`, use a “standard” sample (a Mersenne Twister pseudo-random sample if `flags\$pseudo.random=1`, otherwise a Sobol quasi-random sample). `key2` integer that determines sampling in the final integration phase: same as `key1`, but here \code{n = |key2|} determines the number of points, \code{n > 39}, sample n points, \code{n < 40}, sample \code{n} `nneed` points, where `nneed` is the number of points needed to reach the prescribed accuracy, as estimated by Divonne from the results of the partitioning phase. `key3` integer that sets the strategy for the refinement phase: `key3 = 0`, do not treat the subregion any further. `key3 = 1`, split the subregion up once more. Otherwise, the subregion is sampled a third time with `key3` specifying the sampling parameters exactly as `key2` above. `max.pass` 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 `max.pass` successive iterations. A decrease in points generally indicates that Divonne discovered new structures of the integrand and was able to find a more effective partitioning. `max.pass` can be understood as the number of “safety” iterations that are performed before the partition is accepted as final and counting consequently restarts at zero whenever new structures are found. `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 non-zero `border` if the integrand subroutine cannot produce values directly on the integration boundary. The relative width of the border is identical in all the dimensions. For example, set `border=0.1` for a border of width equal to 10% of the width of the integration region. `max.chisq` the maximum Chi2 value a single subregion is allowed to have in the final integration phase. Regions which fail this Chi2 test and whose sample averages differ by more than `min.deviation` move on to the refinement phase. `min.deviation` 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 Chi2 test. Only if the two sampling averages obtained for the region differ by more than this bound is the region further treated. `xgiven` a matrix ( `ndim`, `ngiven`). A list of `ngiven` points where the integrand might have peaks. Divonne will consider these points when partitioning the integration region. The idea here is to help the integrator find the extrema of the integrand in the presence of very narrow peaks. Even if only the approximate location of such peaks is known, this can considerably speed up convergence. `nextra` the maximum number of extra points the peak-finder subroutine will return. If `nextra` is zero, `peakfinder` is not called and an arbitrary object may be passed in its place, e.g. just 0. `peakfinder` the peak-finder 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 `xgiven`. It is expected to be declared as `peakfinder <- function(bounds)` where `bounds` is a matrix of dimension (`ndim, 2`) which contains the upper and lower bounds of the subregion. The names of the columns are `c("lower", "upper")`. The returned value should be a matrix (`ndim, nx`) where `nx` is the actual number of points (should be less or equal to `nextra`).

## Details

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 half-range).

See details in the documentation.

## Value

Idem as `cuhre`. Here `ifail` may be `>1` when the accuracy goal was not met within the allowed maximum number of integrand evaluations. Divonne can estimate the number of points by which `maxeval` needs to be increased to reach the desired accuracy and returns this value.

## References

J. H. Friedman, M. H. Wright (1981) A nested partitioning procedure for numerical multiple integration. ACM Trans. Math. Software, 7(1), 76-92.

J. H. Friedman, M. H. Wright (1981) User's guide for DIVONNE. SLAC Report CGTM-193-REV, CGTM-193, Stanford University.

T. Hahn (2005) CUBA-a library for multidimensional numerical integration. Computer Physics Communications, 168, 78-95.

`cuhre`, `suave`, `vegas`

## Examples

 ``` 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31``` ```NDIM <- 3 NCOMP <- 1 integrand <- function(arg, phase) { x <- arg[1] y <- arg[2] z <- arg[3] ff <- sin(x)*cos(y)*exp(z); return(ff) } divonne(NDIM, NCOMP, integrand, rel.tol=1e-3, abs.tol=1e-12, flags=list(verbose=2), key1= 47) # Example with a peak-finder function NMAX <- 4 peakf <- function(bounds) { # 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 2:NMAX) { x[,i] <- x[,(i-1)] + pas } return(x) } #end peakf divonne(NDIM, NCOMP, integrand, flags=list(verbose=0) , peakfinder=peakf, nextra=NMAX) ```

### Example output

```Divonne input parameters:
ndim 3
ncomp 1
rel.tol 0.001
abs.tol 1e-12
pseudo.random  0
final 0
verbose 2
min.eval 0
max.eval 50000
key1 47
key2 1
key3 1
max.pass 5
border 0
max.chisq 10
min.deviation 0.25
ngiven 0
nextra 0
Partitioning phase:
Iteration 1 (pass 0):  8 regions
836 integrand evaluations so far,
406 in optimizing regions,
70 in finding cuts
[1] 0.665011 +- 0.00470198
Iteration 2 (pass 0):  9 regions
966 integrand evaluations so far,
478 in optimizing regions,
80 in finding cuts
[1] 0.664964 +- 0.00429467
Iteration 3 (pass 1):  10 regions
1096 integrand evaluations so far,
550 in optimizing regions,
90 in finding cuts
[1] 0.664949 +- 0.00388393
Iteration 4 (pass 2):  11 regions
1194 integrand evaluations so far,
590 in optimizing regions,
100 in finding cuts
[1] 0.664887 +- 0.0035842
Iteration 5 (pass 3):  12 regions
1308 integrand evaluations so far,
646 in optimizing regions,
110 in finding cuts
[1] 0.664853 +- 0.0033385
Iteration 6 (pass 4):  13 regions
1438 integrand evaluations so far,
718 in optimizing regions,
120 in finding cuts
[1] 0.664825 +- 0.0031276
Iteration 7 (pass 5):  14 regions
1568 integrand evaluations so far,
790 in optimizing regions,
130 in finding cuts
[1] 0.664816 +- 0.00292074

Main integration on 14 regions with 211 samples per region.integral: 0.6646195 (+-0.00064)
nregions: 14; number of evaluations:  3052; probability:  1.110223e-16
integral: 0.6646195 (+-0.00064)
nregions: 14; number of evaluations:  3136; probability:  1.110223e-16
```

R2Cuba documentation built on May 29, 2017, 7:53 p.m.