# Version 2.0 Notes

#### Balasubramanian Narasimhan

#### 2025-05-15

## 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, and in some cases, provide somewhat different answers. That still makes sense and depends on the underlying methodology used.

## 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))

    ## Cuhre input parameters:
    ##   ndim 3
    ##   ncomp 1
    ##   nvec 1
    ##   epsrel 1e-05
    ##   epsabs 1e-12
    ##   flags 6
    ##   mineval 0
    ##   maxeval 1000000
    ##   key 0
    ##   statefile "(null)"
    ##
    ## Iteration 1:  127 integrand evaluations so far
    ## [1] 0.330292 +- 0.16199      chisq 0 (0 df)
    ##
    ## Iteration 2:  381 integrand evaluations so far
    ## [1] 0.334551 +- 0.0726708    chisq 0.000575398 (1 df)
    ##
    ## Iteration 3:  635 integrand evaluations so far
    ## [1] 0.334084 +- 0.0126583    chisq 0.000588724 (2 df)
    ##
    ## Iteration 4:  889 integrand evaluations so far
    ## [1] 0.334093 +- 0.00467064   chisq 0.000590497 (3 df)
    ##
    ## Iteration 5:  1143 integrand evaluations so far
    ## [1] 0.33411 +- 0.00213905    chisq 0.00060536 (4 df)
    ##
    ## Iteration 6:  1397 integrand evaluations so far
    ## [1] 0.334113 +- 0.000658593      chisq 0.000616501 (5 df)
    ##
    ## Iteration 7:  1651 integrand evaluations so far
    ## [1] 0.334113 +- 0.000290809      chisq 0.000617561 (6 df)
    ##
    ## Iteration 8:  1905 integrand evaluations so far
    ## [1] 0.334113 +- 4.6551e-05   chisq 0.000617624 (7 df)
    ##
    ## Iteration 9:  2159 integrand evaluations so far
    ## [1] 0.334113 +- 3.69574e-05      chisq 0.000617628 (8 df)
    ##
    ## Iteration 10:  2413 integrand evaluations so far
    ## [1] 0.334113 +- 6.66791e-05      chisq 0.000623038 (9 df)
    ##
    ## Iteration 11:  2667 integrand evaluations so far
    ## [1] 0.334113 +- 2.65953e-05      chisq 0.000637044 (10 df)
    ##
    ## Iteration 12:  2921 integrand evaluations so far
    ## [1] 0.334113 +- 3.88133e-05      chisq 0.000639547 (11 df)
    ##
    ## Iteration 13:  3175 integrand evaluations so far
    ## [1] 0.334113 +- 2.10511e-05      chisq 0.000643648 (12 df)
    ##
    ## Iteration 14:  3429 integrand evaluations so far
    ## [1] 0.334113 +- 1.31307e-05      chisq 0.000650986 (13 df)
    ##
    ## Iteration 15:  3683 integrand evaluations so far
    ## [1] 0.334113 +- 8.69938e-06      chisq 0.000654881 (14 df)
    ##
    ## Iteration 16:  3937 integrand evaluations so far
    ## [1] 0.334113 +- 1.88417e-05      chisq 0.000655658 (15 df)
    ##
    ## Iteration 17:  4191 integrand evaluations so far
    ## [1] 0.334113 +- 1.59906e-05      chisq 0.000656618 (16 df)
    ##
    ## Iteration 18:  4445 integrand evaluations so far
    ## [1] 0.334113 +- 4.85581e-06      chisq 0.000660585 (17 df)
    ##
    ## Iteration 19:  4699 integrand evaluations so far
    ## [1] 0.334113 +- 2.22627e-06      chisq 0.000668304 (18 df)

    ## $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.

    str(default_args())

    ## 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

## 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.5.0 (2025-04-11)
    ## Platform: aarch64-apple-darwin20
    ## Running under: macOS Sequoia 15.4.1
    ##
    ## Matrix products: default
    ## BLAS:   /Library/Frameworks/R.framework/Versions/4.5-arm64/Resources/lib/libRblas.0.dylib
    ## LAPACK: /Library/Frameworks/R.framework/Versions/4.5-arm64/Resources/lib/libRlapack.dylib;  LAPACK version 3.12.1
    ##
    ## 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.3-3  cubature_2.1.3 benchr_0.2.5
    ##
    ## loaded via a namespace (and not attached):
    ##  [1] digest_0.6.37      R6_2.6.1           fastmap_1.2.0      xfun_0.52
    ##  [5] cachem_1.1.0       knitr_1.50         htmltools_0.5.8.1  rmarkdown_2.29
    ##  [9] lifecycle_1.0.4    cli_3.6.5          RcppProgress_0.4.2 sass_0.4.10
    ## [13] jquerylib_0.1.4    compiler_4.5.0     tools_4.5.0        evaluate_1.0.3
    ## [17] bslib_0.9.0        Rcpp_1.0.14.12     yaml_2.3.10        rlang_1.1.6
    ## [21] jsonlite_2.0.0