File: hcubature.Rd

package info (click to toggle)
r-cran-cubature 1.3-6-1
  • links: PTS
  • area: main
  • in suites: stretch
  • size: 6,088 kB
  • ctags: 163
  • sloc: ansic: 1,499; cpp: 81; makefile: 6
file content (353 lines) | stat: -rw-r--r-- 12,022 bytes parent folder | download
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
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
123
124
125
126
127
128
129
130
131
132
133
134
135
136
137
138
139
140
141
142
143
144
145
146
147
148
149
150
151
152
153
154
155
156
157
158
159
160
161
162
163
164
165
166
167
168
169
170
171
172
173
174
175
176
177
178
179
180
181
182
183
184
185
186
187
188
189
190
191
192
193
194
195
196
197
198
199
200
201
202
203
204
205
206
207
208
209
210
211
212
213
214
215
216
217
218
219
220
221
222
223
224
225
226
227
228
229
230
231
232
233
234
235
236
237
238
239
240
241
242
243
244
245
246
247
248
249
250
251
252
253
254
255
256
257
258
259
260
261
262
263
264
265
266
267
268
269
270
271
272
273
274
275
276
277
278
279
280
281
282
283
284
285
286
287
288
289
290
291
292
293
294
295
296
297
298
299
300
301
302
303
304
305
306
307
308
309
310
311
312
313
314
315
316
317
318
319
320
321
322
323
324
325
326
327
328
329
330
331
332
333
334
335
336
337
338
339
340
341
342
343
344
345
346
347
348
349
350
351
352
353
% Generated by roxygen2: do not edit by hand
% Please edit documentation in R/hcubature.R
\name{hcubature}
\alias{adaptIntegrate}
\alias{hcubature}
\alias{pcubature}
\title{Adaptive multivariate integration over hypercubes (hcubature and pcubature)}
\usage{
hcubature(f, lowerLimit, upperLimit, ..., tol = 1e-05, fDim = 1,
  maxEval = 0, absError = 0, doChecking = FALSE,
  vectorInterface = FALSE, norm = c("INDIVIDUAL", "PAIRED", "L2", "L1",
  "LINF"))

pcubature(f, lowerLimit, upperLimit, ..., tol = 1e-05, fDim = 1,
  maxEval = 0, absError = 0, doChecking = FALSE,
  vectorInterface = FALSE, norm = c("INDIVIDUAL", "PAIRED", "L2", "L1",
  "LINF"))
}
\arguments{
\item{f}{The function (integrand) to be integrated}

\item{lowerLimit}{The lower limit of integration, a vector for hypercubes}

\item{upperLimit}{The upper limit of integration, a vector for hypercubes}

\item{...}{All other arguments passed to the function f}

\item{tol}{The maximum tolerance, default 1e-5.}

\item{fDim}{The dimension of the integrand, default 1, bears no relation to
the dimension of the hypercube}

\item{maxEval}{The maximum number of function evaluations needed, default 0
implying no limit. Note that the actual number of function evaluations
performed is only approximately guaranteed not to exceed this number.}

\item{absError}{The maximum absolute error tolerated}

\item{doChecking}{A flag to be a bit anal about checking inputs to C
routines. A FALSE value results in approximately 9 percent speed gain in our
experiments. Your mileage will of course vary. Default value is FALSE.}

\item{vectorInterface}{A flag that indicates whether to use the vector interface
and is by default FALSE. See details below}

\item{norm}{For vector-valued integrands, \code{norm} specifies the norm that is
used to measure the error and determine convergence properties. See below.}
}
\value{
The returned value is a list of three items: \item{integral}{the
value of the integral} \item{error}{the estimated relative error}
\item{functionEvaluations}{the number of times the function was evaluated}
\item{returnCode}{the actual integer return code of the C routine}
}
\description{
The function performs adaptive multidimensional integration (cubature) of
(possibly) vector-valued integrands over hypercubes. The function includes
a vector interface where the integrand may be evaluated at several hundred
points in a single call.
}
\details{
The function merely calls Johnson's C code and returns the results.

One can specify a maximum number of function evaluations (default is 0 for
no limit).  Otherwise, the integration stops when the estimated error is
less than the absolute error requested, or when the estimated error is less
than tol times the integral, in absolute value, or the maximum number of
iterations is reached (see parameter info below), whichever is earlier.

For compatibility with earlier versions, the \code{adaptIntegrate} function
is an alias for the underlying \code{hcubature} function which uses h-adaptive
integration. Otherwise, the calling conventions are the same.

We highly recommend referring to the vignette to achieve the best results!



The \code{hcubature} function is the h-adaptive version that recursively partitions
the integration domain into smaller subdomains, applying the same integration
rule to each, until convergence is achieved.

The p-adaptive version, \code{pcubature}, repeatedly doubles the degree
of the quadrature rules until convergence is achieved, and is based on a tensor
product of Clenshaw-Curtis quadrature rules. This algorithm is often superior
to h-adaptive integration for smooth integrands in a few (<=3) dimensions,
but is a poor choice in higher dimensions or for non-smooth integrands.
Compare with \code{hcubature} which also takes the same arguments.

The vector interface requires the integrand to take a matrix as its argument.
The return value should also be a matrix. The number of points at which the
integrand may be evaluated is not under user control: the integration routine
takes care of that and this number may run to several hundreds. We strongly
advise vectorization; see vignette.

The \code{norm} argument is irrelevant for scalar integrands and is ignored.
Given vectors \eqn{v} and \eqn{e} of estimated integrals and errors therein,
respectively, the \code{norm} argument takes on one of the following values:
\describe{
  \item{\code{INDIVIDUAL}}{Convergence is achieved only when each integrand
  (each component of \eqn{v} and \eqn{e}) individually satisfies the requested
  error tolerances}
  \item{\code{L1}, \code{L2}, \code{LINF}}{The absolute error is measured as
  \eqn{|e|} and the relative error as \eqn{|e|/|v|}, where \eqn{|...|} is the
  \eqn{L_1}, \eqn{L_2}, or \eqn{\L_{\infty}} norm, respectively}
  \item{\code{PAIRED}}{Like \code{INDIVIDUAL}, except that the integrands are
  grouped into consecutive pairs, with the error tolerance applied in an
  \eqn{L_2} sense to each pair. This option is mainly useful for integrating
  vectors of complex numbers, where each consecutive pair of real integrands
  is the real and imaginary parts of a single complex integrand, and the concern
  is only the error in the complex plane rather than the error in the real and
  imaginary parts separately}
}
}
\examples{

\dontrun{
## Test function 0
## Compare with original cubature result of
## ./cubature_test 2 1e-4 0 0
## 2-dim integral, tolerance = 0.0001
## integrand 0: integral = 0.708073, est err = 1.70943e-05, true err = 7.69005e-09
## #evals = 17

testFn0 <- function(x) {
  prod(cos(x))
}

hcubature(testFn0, rep(0,2), rep(1,2), tol=1e-4)

pcubature(testFn0, rep(0,2), rep(1,2), tol=1e-4)

M_2_SQRTPI <- 2/sqrt(pi)

## Test function 1
## Compare with original cubature result of
## ./cubature_test 3 1e-4 1 0
## 3-dim integral, tolerance = 0.0001
## integrand 1: integral = 1.00001, est err = 9.67798e-05, true err = 9.76919e-06
## #evals = 5115

testFn1 <- function(x) {
  val <- sum (((1-x) / x)^2)
  scale <- prod(M_2_SQRTPI/x^2)
  exp(-val) * scale
}

hcubature(testFn1, rep(0, 3), rep(1, 3), tol=1e-4)
pcubature(testFn1, rep(0, 3), rep(1, 3), tol=1e-4)

##
## Test function 2
## Compare with original cubature result of
## ./cubature_test 2 1e-4 2 0
## 2-dim integral, tolerance = 0.0001
## integrand 2: integral = 0.19728, est err = 1.97261e-05, true err = 4.58316e-05
## #evals = 166141

testFn2 <- function(x) {
  ## discontinuous objective: volume of hypersphere
  radius <- as.double(0.50124145262344534123412)
  ifelse(sum(x*x) < radius*radius, 1, 0)
}

hcubature(testFn2, rep(0, 2), rep(1, 2), tol=1e-4)
pcubature(testFn2, rep(0, 2), rep(1, 2), tol=1e-4)

##
## Test function 3
## Compare with original cubature result of
## ./cubature_test 3 1e-4 3 0
## 3-dim integral, tolerance = 0.0001
## integrand 3: integral = 1, est err = 0, true err = 2.22045e-16
## #evals = 33

testFn3 <- function(x) {
  prod(2*x)
}

hcubature(testFn3, rep(0,3), rep(1,3), tol=1e-4)
pcubature(testFn3, rep(0,3), rep(1,3), tol=1e-4)

##
## Test function 4 (Gaussian centered at 1/2)
## Compare with original cubature result of
## ./cubature_test 2 1e-4 4 0
## 2-dim integral, tolerance = 0.0001
## integrand 4: integral = 1, est err = 9.84399e-05, true err = 2.78894e-06
## #evals = 1853

testFn4 <- function(x) {
  a <- 0.1
  s <- sum((x - 0.5)^2)
  (M_2_SQRTPI / (2. * a))^length(x) * exp (-s / (a * a))
}

hcubature(testFn4, rep(0,2), rep(1,2), tol=1e-4)
pcubature(testFn4, rep(0,2), rep(1,2), tol=1e-4)

##
## Test function 5 (double Gaussian)
## Compare with original cubature result of
## ./cubature_test 3 1e-4 5 0
## 3-dim integral, tolerance = 0.0001
## integrand 5: integral = 0.999994, est err = 9.98015e-05, true err = 6.33407e-06
## #evals = 59631

testFn5 <- function(x) {
  a <- 0.1
  s1 <- sum((x - 1/3)^2)
  s2 <- sum((x - 2/3)^2)
  0.5 * (M_2_SQRTPI / (2. * a))^length(x) * (exp(-s1 / (a * a)) + exp(-s2 / (a * a)))
}

hcubature(testFn5, rep(0,3), rep(1,3), tol=1e-4)
pcubature(testFn5, rep(0,3), rep(1,3), tol=1e-4)

##
## Test function 6 (Tsuda's example)
## Compare with original cubature result of
## ./cubature_test 4 1e-4 6 0
## 4-dim integral, tolerance = 0.0001
## integrand 6: integral = 0.999998, est err = 9.99685e-05, true err = 1.5717e-06
## #evals = 18753

testFn6 <- function(x) {
  a <- (1 + sqrt(10.0)) / 9.0
  prod(a / (a + 1) * ((a + 1) / (a + x))^2)
}

hcubature(testFn6, rep(0,4), rep(1,4), tol=1e-4)
pcubature(testFn6, rep(0,4), rep(1,4), tol=1e-4)


##
## Test function 7
##   test integrand from W. J. Morokoff and R. E. Caflisch, "Quasi=
##   Monte Carlo integration," J. Comput. Phys 122, 218-230 (1995).
##   Designed for integration on [0,1]^dim, integral = 1. */
## Compare with original cubature result of
## ./cubature_test 3 1e-4 7 0
## 3-dim integral, tolerance = 0.0001
## integrand 7: integral = 1.00001, est err = 9.96657e-05, true err = 1.15994e-05
## #evals = 7887

testFn7 <- function(x) {
  n <- length(x)
  p <- 1/n
  (1 + p)^n * prod(x^p)
}

hcubature(testFn7, rep(0,3), rep(1,3), tol=1e-4)
pcubature(testFn7, rep(0,3), rep(1,3), tol=1e-4)


## Example from web page
## http://ab-initio.mit.edu/wiki/index.php/Cubature
##
## f(x) = exp(-0.5(euclidean_norm(x)^2)) over the three-dimensional
## hyperbcube [-2, 2]^3
## Compare with original cubature result
testFnWeb <-  function(x) {
  exp(-0.5 * sum(x^2))
}

hcubature(testFnWeb, rep(-2,3), rep(2,3), tol=1e-4)
pcubature(testFnWeb, rep(-2,3), rep(2,3), tol=1e-4)

## Test function I.1d from
## Numerical integration using Wang-Landau sampling
## Y. W. Li, T. Wust, D. P. Landau, H. Q. Lin
## Computer Physics Communications, 2007, 524-529
## Compare with exact answer: 1.63564436296
##
I.1d <- function(x) {
  sin(4*x) *
    x * ((x * ( x * (x*x-4) + 1) - 1))
}

hcubature(I.1d, -2, 2, tol=1e-7)
pcubature(I.1d, -2, 2, tol=1e-7)

## Test function I.2d from
## Numerical integration using Wang-Landau sampling
## Y. W. Li, T. Wust, D. P. Landau, H. Q. Lin
## Computer Physics Communications, 2007, 524-529
## Compare with exact answer: -0.01797992646
##
##
I.2d <- function(x) {
  x1 = x[1]
  x2 = x[2]
  sin(4*x1+1) * cos(4*x2) * x1 * (x1*(x1*x1)^2 - x2*(x2*x2 - x1) +2)
}

hcubature(I.2d, rep(-1, 2), rep(1, 2), maxEval=10000)
pcubature(I.2d, rep(-1, 2), rep(1, 2), maxEval=10000)

##
## Example of multivariate normal integration borrowed from
## package mvtnorm (on CRAN) to check ... argument
## Compare with output of
## pmvnorm(lower=rep(-0.5, m), upper=c(1,4,2), mean=rep(0, m), corr=sigma, alg=Miwa())
##     0.3341125.  Blazing quick as well!  Ours is, not unexpectedly, much slower.
##
dmvnorm <- function (x, mean, sigma, log = FALSE) {
    if (is.vector(x)) {
        x <- matrix(x, ncol = length(x))
    }
    if (missing(mean)) {
        mean <- rep(0, length = ncol(x))
    }
    if (missing(sigma)) {
        sigma <- diag(ncol(x))
    }
    if (NCOL(x) != NCOL(sigma)) {
        stop("x and sigma have non-conforming size")
    }
    if (!isSymmetric(sigma, tol = sqrt(.Machine$double.eps),
        check.attributes = FALSE)) {
        stop("sigma must be a symmetric matrix")
    }
    if (length(mean) != NROW(sigma)) {
        stop("mean and sigma have non-conforming size")
    }
    distval <- mahalanobis(x, center = mean, cov = sigma)
    logdet <- sum(log(eigen(sigma, symmetric = TRUE, only.values = TRUE)$values))
    logretval <- -(ncol(x) * log(2 * pi) + logdet + distval)/2
    if (log)
        return(logretval)
    exp(logretval)
}

m <- 3
sigma <- diag(3)
sigma[2,1] <- sigma[1, 2] <- 3/5 ; sigma[3,1] <- sigma[1, 3] <- 1/3
sigma[3,2] <- sigma[2, 3] <- 11/15
hcubature(dmvnorm, lower=rep(-0.5, m), upper=c(1,4,2),
                        mean=rep(0, m), sigma=sigma, log=FALSE,
               maxEval=10000)
pcubature(dmvnorm, lower=rep(-0.5, m), upper=c(1,4,2),
                        mean=rep(0, m), sigma=sigma, log=FALSE,
               maxEval=10000)
}

}
\author{
Balasubramanian Narasimhan
}
\references{
See \url{http://ab-initio.mit.edu/wiki/index.php/Cubature}.
}
\keyword{math}