Inaccurate approximation and some ideas

[jobstdavid]

Hey,

I tested your R-package "approxcdf" and I like it a lot, especially due to its speed! That's great! Unfortunately, I detected a completely inaccurate approximation for the trivariate normal distribution:

library(approxcdf)
library(mvtnorm)
set.seed(435)
Cov <- randcorr::randcorr(3)
upper <- c(0, 0, Inf)
# result 0
approxcdf::pmvn(upper, Cov = Cov, mean = rep(0, 3))
# result 0.1551707
as.numeric(mvtnorm::pmvnorm(upper = upper, corr = Cov, algorithm = "TVPACK"))
as.numeric(mvtnorm::pmvnorm(upper = upper, corr = Cov, algorithm = "GenzBretz"))

# hacky solution (result 0.1551707)
approxcdf::pmvn(c(0, 0, 10000), Cov = Cov, mean = rep(0, 3))

This issue also happens for

upper <- c(0, Inf, Inf)

i.e. settings where more than one upper bound is infinity are affected, as well. Do you think, you can fix it?

Some other ideas came to my mind, which could improve the ease of use and you could think about it:

• Separate functions for the bivariate and trivariate normal distribution: The two-dimensional and three-dimensional case are very often of interest for the user.
• Vectorized functions with arguments (x1, x2, rho) or (x1, x2, x3, rho12, rho13, rho23) as called in https://github.com/david-cortes/approxcdf/blob/ff51ad1a177c4be7fc4902f60c6947e2ae3da5b0/src/genz.cpp might provide a further advantage over mvtnorm and relates to the first point.
• Have you thought about publishing the package on CRAN?

[david-cortes]

Thanks for pointing this out.

• The methods as implemented in this package do not work with infinite values. I've added a note in the docs and extra check.
• While there isn't any dedicated ptvn function offered here, it does have special cases for handling n<=3: https://github.com/david-cortes/approxcdf/blob/ff51ad1a177c4be7fc4902f60c6947e2ae3da5b0/src/tvbs.cpp#L198 The method for 3d is very similar to TVPACK. Note that while it is theoretically possible to vectorize and it has compiler hints for it, it's very unlikely that a compiler will actually end up vectorizing that (at least GCC and CLANG do not), since it includes operations like exp and sin.
• I don't feel this type of package belongs in CRAN - there's plenty of packages offering more exact CDF calculations with varying degrees of speed, and approximate calculations as implemented here could be very harmful to use in place of exact calculations in other packages. I wouldn't want to find out that some model from a package I'm using produces incorrect results due to inaccurate CDF calculations for example. For more serious mass usage, chances are you'd want to use the C interface to integrate it somewhere rather than the R/Python package.

[david-cortes]

Also, if you're looking for a fast bivariate/trivariate CDF calculator to use in some package, you can shave off most of the running time latency from mvtnorm by calling the TVPACK Fortran functions directly without the extra checks that mvtnorm does: https://www.math.wsu.edu/faculty/genz/software/software.html

Timings from it are not too much worse than from the implementations in this package, but the calculations are more accurate.

[jobstdavid]

@david-cortes Thank you for your helpful hints!