My goal is to find a faster way to calculate something like
mvtnorm::pmvnorm(upper = rep(1,100))
that is, the tail probability of multivariate normal distribution with mean 0 and arbitrary covariance matrix. The
upper bound is also arbitrary.
pmvnorm implements three algorithms:
GenzBretz (up to dimension 1000),
Miwa (up to dimension 20) and
TVPACK (2- and 3-dimensional problems for semi-infinite integration regions). For my case of dimension 100, I can only use
GenzBretz and it is the golden standard. However, for each run of the above code, it would take 0.6 second and it is called 40 times in each of 600 iterations. Thus the whole algorithm takes a long time.
To speed up the computation, I have tried the following strategies:
Decrease the accuracy of
GenzBretz. I think
abseps(absolute error tolerance) is the key parameter, but the default seems already liberal and relaxing it does not improve the speed largely.
GenzBretz(maxpts = 25000, abseps = 0.001, releps = 0)
Find alternative functions. I found
OpenMx::omxMnoris faster than
pmvnormin low dimensions, but it is much slower when the dimension becomes as high as 100. Another function is
bayesm::ghkvec, which computes the so called GHK approximation using an importance sampling method. So the tricky part is that the result is not fixed, which causes problems in my algorithm.
Search the literature. I tested the approach described in this paper but found it is slower than
Develop a customized approximation. I was suggested to use Laplace approximation and quadrature methods. But I did not find a readily available function for the multivariate case with semi-infinite bounds. It is also not promising in theory since there is no N data points.
So my question is, is there any function faster than
pmvnorm when the dimension is 100, and can maintain a (fixed and) relatively accurate approximation for multivariate normal probabilities?