# Normal multivariate orthant probabilities

(Previously I posted a similar question on math.SE, hoping that this question would have an easy answer. As the question appears hard, I am hoping I can perhaps get more feedback here.)

Let $$\mathbf{X} \sim N(\mathbf{0}, \mathbf{\Sigma})$$ be a $$k$$-dimensional Gaussian vector with non-trivial covariance matrix $$\mathbf{\Sigma}$$. I am interested in $$p = \Pr(\mathbf{X} \geq \mathbf{0})$$, i.e. the probability that all coordinates are non-negative. This is also known as the orthant probability for $$\mathbf{X}$$, and explicit formulas for orthant probabilities for arbitrary covariance matrices $$\mathbf{\Sigma}$$ are known in $$2$$ and $$3$$ dimensions, with some special cases having been studied in $$4$$ dimensions, and very few explicit formulas apparently known for higher-dimensional cases. I am interested in these orthant probabilities for various different matrices $$\mathbf{\Sigma}$$, perhaps the simplest of which is the following: $$\mathbf{\Sigma} = \begin{pmatrix} 1 & u \\ u & 1 \end{pmatrix} \otimes \begin{pmatrix} 1 & 1/2 & 1/2 \\ 1/2 & 1 & 1/2 \\ 1/2 & 1/2 & 1 \end{pmatrix} = \begin{pmatrix} 1 & 1/2 & 1/2 & u & u/2 & u/2 \\ 1/2 & 1 & 1/2 & u/2 & u & u/2 \\ 1/2 & 1/2 & 1 & u/2 & u/2 & u \\ u & u/2 & u/2 & 1 & 1/2 & 1/2 \\ u/2 & u & u/2 & 1/2 & 1 & 1/2 \\ u/2 & u/2 & u & 1/2 & 1/2 & 1 \end{pmatrix}$$ Here $$u$$ is some constant between $$0$$ and $$1$$.

Most of the literature I found on this topic (for actually deriving closed-form expressions for these probabilities) dates back to long ago, e.g. works by Cheng, Childs, David, Plackett, Steck in the 1950s and 1960s (mostly in the journal Biometrika). I am not so hopeful, but it would be great if

• someone could find a closed form for the above case;
• someone could point me to literature I might have missed on finding closed form expressions for such high-dimensional cases;
• someone could explain which techniques/strategies may generally be useful for finding closed-form expressions for even dimensions.

As for a slightly fishier approach, which nonetheless might be rewarding: I know that orthant probabilities for a few low-dimensional cases can be written as follows: \begin{align} \mathbf{\Sigma} = \begin{pmatrix} 1 & a \\ a & 1 \end{pmatrix} &\implies p = \frac{\arccos (-a)}{2 \pi}, \\ \mathbf{\Sigma} = \begin{pmatrix} 1 & a & b \\ a & 1 & c \\ b & c & 1 \end{pmatrix} &\implies p = \frac{\arccos (-a)}{4 \pi} + \frac{\arccos (-c)}{4 \pi} - \frac{\arccos (b)}{4 \pi}, \\ \mathbf{\Sigma} = \begin{pmatrix} 1 & a & b & ab \\ a & 1 & ab & b \\ b & ab & 1 & a \\ ab & b & a & 1 \end{pmatrix} &\implies p = \left(\frac{\arccos (-a)}{2 \pi}\right)^2 + \left(\frac{\arccos (-b)}{2 \pi}\right)^2 - \left(\frac{\arccos (a b)}{2 \pi}\right)^2. \end{align} Note that the latter case corresponds to: \begin{align} \mathbf{\Sigma} = \begin{pmatrix} 1 & a \\ a & 1 \end{pmatrix} \otimes \begin{pmatrix} 1 & b \\ b & 1 \end{pmatrix} = \begin{pmatrix} 1 & a & b & ab \\ a & 1 & ab & b \\ b & ab & 1 & a \\ ab & b & a & 1 \end{pmatrix} \end{align} In these cases at least, the probability can be written as a combination of arccosines of the off-diagonal entries, sometimes with minus signs. Alternatively, for the cases of two and four dimensions, we can equivalently express the probabilities in the off-diagonal entries (without changing signs) by considering the inverse matrix $$\mathbf{\Sigma}^{-1}$$ instead: \begin{align} \mathbf{\Sigma}^{-1} = \begin{pmatrix} 1 & -a \\ -a & 1 \end{pmatrix} &\implies p = \frac{\arccos (-a)}{2 \pi}, \\ \mathbf{\Sigma}^{-1} = \begin{pmatrix} 1 & -a & -b & ab \\ -a & 1 & ab & -b \\ -b & ab & 1 & -a \\ ab & -b & -a & 1 \end{pmatrix} &\implies p = \left(\frac{\arccos (-a)}{2 \pi}\right)^2 + \left(\frac{\arccos (-b)}{2 \pi}\right)^2 - \left(\frac{\arccos (a b)}{2 \pi}\right)^2. \end{align} Can someone perhaps establish/conjecture a pattern that might allow to extrapolate to higher-dimensional cases?

Edit 1: I fixed a mistake in the formula for three dimensions. For references: the two- and three-dimensional cases can, among others, be found in the one-page paper Dav53 or the paper Chi67, while the four-dimensional case stated above appears in the appendix of Che68.

Edit 2: For the six-dimensional case, substituting $$u = 1/2$$ leads to $$p \approx 0.115024$$. The six decimals seem quite precise; independently, both $$10^8$$ Monte Carlo experiments in C, and using numerical integration in R, I get these six decimals. So the function $$p(u)$$ should satisfy $$p(0) = 0.25$$, $$p(1) = 0.0625$$, and $$p(1/2) \approx 0.115024$$.

• @MattF. Are you sure it is not $1/4$ and $1/16$? (all positive = 1/2, all negative =(flip) all positive =1/2; what chance do all other cases have after that?) – fedja Jun 17 '19 at 21:55
• @MattF. Not sure about the formulae, but if $c=1$, we just have $3$ Gaussian r.v. with a non-degenerate covariance matrix, so the whole space is realizable while you claim that the Gaussian distribution lives just in two octants. Doesn't this bother you a little bit? – fedja Jun 17 '19 at 22:08
• @TMM In other words, the original expression was $1/4$ above the truth. OK, now it makes sense. :) – fedja Jun 18 '19 at 0:36
• Note that there are many ways to write the formulas; we can also write it just in terms of $\arccos a, \arccos b, \arccos c$ using $\arccos (-a) = \pi - \arccos a$. For the four-dimensional case though, this seemed to me the "cleanest", as otherwise we would get many non-vanishing cross-terms due to the squarings. – TMM Jun 18 '19 at 0:38