Qudrivariate and pentavariate normal orthant probabilities

Two functions are included. The first is 'quadriorth.m' which calculates the orthants probabilities of the quadrivariate normal distribution. Imagine a 4-dimensional space divided into sub-spaces according to the sign of each variable, which results in 2^4 = 16 orthants, say + + + +, + + + –, + + – +, + – + +, ..., and – – – –. This function gives the probability content in each of these 16 orthants. The method is based on a recent paper by Sinn and Keller (2010) which reduces the 4-dimensional integral to the sum of four one-dimensional ones (which can be treated numerically as one integral).

Similarly, the function 'pentaorth.m' calculates the orthant probabilities of the pentavariate normal distribution. In this case we have 2^5 = 32 orthants, say + + + + +, + + + + –, + + + – +, ..., etc. The calculations are based on a simple formula due to Kendall and Stuart (1976) which express the pentavariate orthant probabilities in terms of lower dimension orthants (quadrivariate, trivariate, bivariate, and univariate). Note that 'pentaorth.m' issues multiple calls to 'quadriorth.m'.

The orthant probabilities in both cases above can be calculated by using the multivariate normal cumulative distribution function 'mvncdf' in the Matlab statistics toolbox. However, the provided two functions are much faster, especially if high accuracy of the results is considered. The code may benefit from some vectorization and consideration of symmetries, which will hopefully be achieved in upcoming versions to further enhance speed.