Version 2.0 integrates two well-known cubature libraries in one place:
It also provides a single function cubintegrate
that
allows one to call all methods in a uniform fashion, as I explain
below.
N.B. One has to be aware that there are cases where one library will integrate a function while the other won’t, and in some cases, provide somewhat different answers. That still makes sense and depends on the underlying methodology used.
Following a suggestion by Simen Guare, we now have a function
cubintegrate
that can be used to try out various
integration methods easily. Some examples.
library(cubature)
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
logdet <- sum(log(eigen(sigma, symmetric = TRUE, only.values = TRUE)$values))
my_dmvnorm <- function (x, mean, sigma, logdet) {
x <- matrix(x, ncol = length(x))
distval <- stats::mahalanobis(x, center = mean, cov = sigma)
exp(-(3 * log(2 * pi) + logdet + distval)/2)
}
First we try the scalar invocation with hcubature
.
cubintegrate(f = my_dmvnorm, lower = rep(-0.5, 3), upper = c(1, 4, 2), method = "pcubature",
mean = rep(0, m), sigma = sigma, logdet = logdet)
## $integral
## [1] 0.3341125
##
## $error
## [1] 2.21207e-06
##
## $neval
## [1] 4913
##
## $returnCode
## [1] 0
We can compare that with Cuba’s cuhre
.
cubintegrate(f = my_dmvnorm, lower = rep(-0.5, 3), upper = c(1, 4, 2), method = "cuhre",
mean = rep(0, m), sigma = sigma, logdet = logdet)
## $integral
## [1] 0.3341125
##
## $error
## [1] 2.226266e-06
##
## $nregions
## [1] 19
##
## $neval
## [1] 4699
##
## $prob
## [1] 0
##
## $returnCode
## [1] 0
The Cuba routine can take various further arguments; see for example,
the help on cuhre
. Such arguments can be directly passed to
cubintegrate
.
cubintegrate(f = my_dmvnorm, lower = rep(-0.5, 3), upper = c(1, 4, 2), method = "cuhre",
mean = rep(0, m), sigma = sigma, logdet = logdet,
flags = list(verbose = 2))
## $integral
## [1] 0.3341125
##
## $error
## [1] 2.226266e-06
##
## $nregions
## [1] 19
##
## $neval
## [1] 4699
##
## $prob
## [1] 0
##
## $returnCode
## [1] 0
As there are many such method-specific arguments, you may find the
function default_args()
useful.
## List of 6
## $ hcubature:List of 1
## ..$ norm: chr [1:5] "INDIVIDUAL" "PAIRED" "L2" "L1" ...
## $ pcubature:List of 1
## ..$ norm: chr [1:5] "INDIVIDUAL" "PAIRED" "L2" "L1" ...
## $ cuhre :List of 4
## ..$ minEval : int 0
## ..$ stateFile: NULL
## ..$ flags :List of 6
## .. ..$ verbose : int 0
## .. ..$ final : int 1
## .. ..$ smooth : int 0
## .. ..$ keep_state: int 0
## .. ..$ load_state: int 0
## .. ..$ level : int 0
## ..$ key : int 0
## $ divonne :List of 14
## ..$ minEval : int 0
## ..$ stateFile : NULL
## ..$ flags :List of 6
## .. ..$ verbose : int 0
## .. ..$ final : int 1
## .. ..$ smooth : int 0
## .. ..$ keep_state: int 0
## .. ..$ load_state: int 0
## .. ..$ level : int 0
## ..$ rngSeed : int 0
## ..$ key1 : int 47
## ..$ key2 : int 1
## ..$ key3 : int 1
## ..$ maxPass : int 5
## ..$ border : num 0
## ..$ maxChisq : num 10
## ..$ minDeviation: num 0.25
## ..$ xGiven : NULL
## ..$ nExtra : int 0
## ..$ peakFinder : NULL
## $ sauve :List of 7
## ..$ minEval : int 0
## ..$ stateFile: NULL
## ..$ flags :List of 6
## .. ..$ verbose : int 0
## .. ..$ final : int 1
## .. ..$ smooth : int 0
## .. ..$ keep_state: int 0
## .. ..$ load_state: int 0
## .. ..$ level : int 0
## ..$ rngSeed : int 0
## ..$ nNew : int 1000
## ..$ nMin : int 50
## ..$ flatness : num 50
## $ vegas :List of 8
## ..$ minEval : int 0
## ..$ stateFile: NULL
## ..$ flags :List of 6
## .. ..$ verbose : int 0
## .. ..$ final : int 1
## .. ..$ smooth : int 0
## .. ..$ keep_state: int 0
## .. ..$ load_state: int 0
## .. ..$ level : int 0
## ..$ rngSeed : int 0
## ..$ nStart : int 1000
## ..$ nIncrease: int 500
## ..$ nBatch : int 1000
## ..$ gridNo : int 0
cubintegrate
provides vector intefaces too: the
parameter nVec
is by default 1, indicating a scalar
interface. Any value > 1 results in a vectorized call. So
f
has to be constructed appropriately, thus:
my_dmvnorm_v <- function (x, mean, sigma, logdet) {
distval <- stats::mahalanobis(t(x), center = mean, cov = sigma)
exp(matrix(-(3 * log(2 * pi) + logdet + distval)/2, ncol = ncol(x)))
}
Here, the two underlying C libraries differ. The cubature library
manages the number of points used in vectorization dynamically and this
number can even vary from call to call. So any value of
nVec
greater than 1 is merely a flag to use vectorization.
The Cuba C library on the other hand, will use the actual value of
nVec
.
cubintegrate(f = my_dmvnorm_v, lower = rep(-0.5, 3), upper = c(1, 4, 2), method = "pcubature",
mean = rep(0, m), sigma = sigma, logdet = logdet,
nVec = 128)
## $integral
## [1] 0.3341125
##
## $error
## [1] 2.21207e-06
##
## $neval
## [1] 4913
##
## $returnCode
## [1] 0
cubintegrate(f = my_dmvnorm_v, lower = rep(-0.5, 3), upper = c(1, 4, 2), method = "cuhre",
mean = rep(0, m), sigma = sigma, logdet = logdet,
nVec = 128)
## $integral
## [1] 0.3341125
##
## $error
## [1] 2.226266e-06
##
## $nregions
## [1] 19
##
## $neval
## [1] 4699
##
## $prob
## [1] 0
##
## $returnCode
## [1] 0
## R version 4.4.2 (2024-10-31)
## Platform: x86_64-pc-linux-gnu
## Running under: Ubuntu 24.04.1 LTS
##
## Matrix products: default
## BLAS: /usr/lib/x86_64-linux-gnu/openblas-pthread/libblas.so.3
## LAPACK: /usr/lib/x86_64-linux-gnu/openblas-pthread/libopenblasp-r0.3.26.so; LAPACK version 3.12.0
##
## locale:
## [1] LC_CTYPE=en_US.UTF-8 LC_NUMERIC=C
## [3] LC_TIME=en_US.UTF-8 LC_COLLATE=C
## [5] LC_MONETARY=en_US.UTF-8 LC_MESSAGES=en_US.UTF-8
## [7] LC_PAPER=en_US.UTF-8 LC_NAME=C
## [9] LC_ADDRESS=C LC_TELEPHONE=C
## [11] LC_MEASUREMENT=en_US.UTF-8 LC_IDENTIFICATION=C
##
## time zone: Etc/UTC
## tzcode source: system (glibc)
##
## attached base packages:
## [1] stats graphics grDevices utils datasets methods base
##
## other attached packages:
## [1] mvtnorm_1.3-2 cubature_2.1.1 benchr_0.2.5
##
## loaded via a namespace (and not attached):
## [1] digest_0.6.37 R6_2.5.1 fastmap_1.2.0 xfun_0.49
## [5] maketools_1.3.1 cachem_1.1.0 knitr_1.49 htmltools_0.5.8.1
## [9] rmarkdown_2.29 buildtools_1.0.0 lifecycle_1.0.4 cli_3.6.3
## [13] RcppProgress_0.4.2 sass_0.4.9 jquerylib_0.1.4 compiler_4.4.2
## [17] sys_3.4.3 tools_4.4.2 evaluate_1.0.1 bslib_0.8.0
## [21] Rcpp_1.0.13-1 yaml_2.3.10 jsonlite_1.8.9 rlang_1.1.4