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.
See Also
cuhre, suave, vegas
Examples
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)
Results
R version 3.3.1 (2016-06-21) -- "Bug in Your Hair"
Copyright (C) 2016 The R Foundation for Statistical Computing
Platform: x86_64-pc-linux-gnu (64-bit)
R is free software and comes with ABSOLUTELY NO WARRANTY.
You are welcome to redistribute it under certain conditions.
Type 'license()' or 'licence()' for distribution details.
R is a collaborative project with many contributors.
Type 'contributors()' for more information and
'citation()' on how to cite R or R packages in publications.
Type 'demo()' for some demos, 'help()' for on-line help, or
'help.start()' for an HTML browser interface to help.
Type 'q()' to quit R.
> library(R2Cuba)
> png(filename="/home/ddbj/snapshot/RGM3/R_CC/result/R2Cuba/divonne.Rd_%03d_medium.png", width=480, height=480)
> ### Name: divonne
> ### Title: Integration by Stratified Sampling for Variance Reduction
> ### Aliases: divonne
> ### Keywords: math
>
> ### ** Examples
>
> 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)
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
>
> # 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)
integral: 0.6646195 (+-0.00064)
nregions: 14; number of evaluations: 3136; probability: 1.110223e-16
>
>
>
>
>
> dev.off()
null device
1
>