Version 2.0 Notes

The hardware and bandwidth for this mirror is donated by dogado GmbH, the Webhosting and Full Service-Cloud Provider. Check out our Wordpress Tutorial.
If you wish to report a bug, or if you are interested in having us mirror your free-software or open-source project, please feel free to contact us at mirror[@]dogado.de.

What’s new

Version 2.0 integrates two well-known cubature libraries in one place:

The cubature C library of Steven G. Johnson.
The Cuba C library of Thomas Hahn.

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. Where possible, I have brought such problems to the notice of the authors underlying C library and I encourage you to do the same, by providing reproducible examples.

Unified Interface

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.

default_args()

## $hcubature
## $hcubature$norm
## [1] "INDIVIDUAL" "PAIRED"     "L2"         "L1"         "LINF"      
## 
## 
## $pcubature
## $pcubature$norm
## [1] "INDIVIDUAL" "PAIRED"     "L2"         "L1"         "LINF"      
## 
## 
## $cuhre
## $cuhre$minEval
## [1] 0
## 
## $cuhre$stateFile
## NULL
## 
## $cuhre$flags
## $cuhre$flags$verbose
## [1] 0
## 
## $cuhre$flags$final
## [1] 1
## 
## $cuhre$flags$smooth
## [1] 0
## 
## $cuhre$flags$keep_state
## [1] 0
## 
## $cuhre$flags$load_state
## [1] 0
## 
## $cuhre$flags$level
## [1] 0
## 
## 
## $cuhre$key
## [1] 0
## 
## 
## $divonne
## $divonne$minEval
## [1] 0
## 
## $divonne$stateFile
## NULL
## 
## $divonne$flags
## $divonne$flags$verbose
## [1] 0
## 
## $divonne$flags$final
## [1] 1
## 
## $divonne$flags$smooth
## [1] 0
## 
## $divonne$flags$keep_state
## [1] 0
## 
## $divonne$flags$load_state
## [1] 0
## 
## $divonne$flags$level
## [1] 0
## 
## 
## $divonne$rngSeed
## [1] 0
## 
## $divonne$key1
## [1] 47
## 
## $divonne$key2
## [1] 1
## 
## $divonne$key3
## [1] 1
## 
## $divonne$maxPass
## [1] 5
## 
## $divonne$border
## [1] 0
## 
## $divonne$maxChisq
## [1] 10
## 
## $divonne$minDeviation
## [1] 0.25
## 
## $divonne$xGiven
## NULL
## 
## $divonne$nExtra
## [1] 0
## 
## $divonne$peakFinder
## NULL
## 
## 
## $sauve
## $sauve$minEval
## [1] 0
## 
## $sauve$stateFile
## NULL
## 
## $sauve$flags
## $sauve$flags$verbose
## [1] 0
## 
## $sauve$flags$final
## [1] 1
## 
## $sauve$flags$smooth
## [1] 0
## 
## $sauve$flags$keep_state
## [1] 0
## 
## $sauve$flags$load_state
## [1] 0
## 
## $sauve$flags$level
## [1] 0
## 
## 
## $sauve$rngSeed
## [1] 0
## 
## $sauve$nNew
## [1] 1000
## 
## $sauve$nMin
## [1] 50
## 
## $sauve$flatness
## [1] 50
## 
## 
## $vegas
## $vegas$minEval
## [1] 0
## 
## $vegas$stateFile
## NULL
## 
## $vegas$flags
## $vegas$flags$verbose
## [1] 0
## 
## $vegas$flags$final
## [1] 1
## 
## $vegas$flags$smooth
## [1] 0
## 
## $vegas$flags$keep_state
## [1] 0
## 
## $vegas$flags$load_state
## [1] 0
## 
## $vegas$flags$level
## [1] 0
## 
## 
## $vegas$rngSeed
## [1] 0
## 
## $vegas$nStart
## [1] 1000
## 
## $vegas$nIncrease
## [1] 500
## 
## $vegas$nBatch
## [1] 1000
## 
## $vegas$gridNo
## [1] 0

Vectorization

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

Session Info

sessionInfo()

## R version 4.3.0 (2023-04-21)
## Platform: x86_64-apple-darwin22.4.0 (64-bit)
## Running under: macOS Ventura 13.4
## 
## Matrix products: default
## BLAS:   /usr/local/Cellar/openblas/0.3.23/lib/libopenblasp-r0.3.23.dylib 
## LAPACK: /usr/local/Cellar/r/4.3.0_1/lib/R/lib/libRlapack.dylib;  LAPACK version 3.11.0
## 
## locale:
## [1] en_US.UTF-8/en_US.UTF-8/en_US.UTF-8/C/en_US.UTF-8/en_US.UTF-8
## 
## time zone: America/Los_Angeles
## tzcode source: internal
## 
## attached base packages:
## [1] stats     graphics  grDevices utils     datasets  methods   base     
## 
## other attached packages:
## [1] mvtnorm_1.1-3  cubature_2.1.0 benchr_0.2.5  
## 
## loaded via a namespace (and not attached):
##  [1] digest_0.6.31      R6_2.5.1           fastmap_1.1.1      xfun_0.39         
##  [5] cachem_1.0.8       knitr_1.43         htmltools_0.5.5    rmarkdown_2.22    
##  [9] cli_3.6.1          RcppProgress_0.4.2 sass_0.4.6         jquerylib_0.1.4   
## [13] compiler_4.3.0     tools_4.3.0        evaluate_0.21      bslib_0.5.0       
## [17] Rcpp_1.0.10        yaml_2.3.7         rlang_1.1.1        jsonlite_1.8.5

These binaries (installable software) and packages are in development.
They may not be fully stable and should be used with caution. We make no claims about them.
Health stats visible at Monitor.

Version 2.0 Notes

Balasubramanian Narasimhan

2023-06-11

What’s new

Unified Interface

Vectorization

Session Info