Normal multivariate orthant probabilities

Question

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

Answer 1

Define higher arcsines as follows:

$$\arcsin(a) = \int_0^a \frac 1 { \sqrt{1-\alpha^2} } d\alpha$$

$$\arcsin(a,b,c) = \int_0^b \frac { \arcsin\left( a \beta c / v / w \right) } { \sqrt{1-\beta^2} } d\beta$$

$$\text{where}\; v = \sqrt{ 1-a^2-\beta^2 },\; w = \sqrt{ 1-\beta^2-c^2 }.$$

$$\arcsin(a,b,c,d,e) = \int_0^c \frac { \arcsin( x )\arcsin( y )+\arcsin( x, b \gamma d / v / w, y ) } { \sqrt{1-\gamma^2} } d\gamma$$

$$\text{where}\; v = \sqrt{1-b^2-\gamma^2},\; w = \sqrt{1-\gamma^2-d^2},\; x = \frac{a}{v} \sqrt{1-\gamma^2},\; y = \frac{e}{w} \sqrt{1-\gamma^2}.$$

Now define the orthoscheme volume:

$$V( a, b, c, d, e ) = \frac {\pi^3} {384} + \frac {\pi^2} {192} \left( \arcsin( a ) + \arcsin( b ) + \arcsin( c ) + \arcsin( d ) + \arcsin( e ) \right) + \frac \pi { 96 } \left( \arcsin( a ) \arcsin( c ) + \arcsin( a ) \arcsin( d ) + \arcsin( a ) \arcsin( e ) + \arcsin( b ) \arcsin( d ) + \arcsin( b ) \arcsin( e ) + \arcsin( c ) \arcsin( e ) + \arcsin( a, b, c ) + \arcsin( b, c, d ) + \arcsin( c, d, e ) \right) + \frac 1 {48} \arcsin( a, b, c, d, e ) .$$

Set the following in terms of parameter $u$:

$$\begin{array}{lll} a = 1-u^2 & b = 1+u^2 & d = 1+2 u^2 \\ e = 3-2 u^2 & f = 3+5 u^2 & g = 1+7 u^2 \\ h = 7-u^2 & i = 3+18 u^2-5 u^4 & j = 2-u^2 \\ k = 11-3 u^2 & l = 3+u^2 & m = 5-u^2 \\ n = 4-3 u^2 & p = 3-u^2 & q = 1+5 u^2 \end{array}.$$

Then for $0 < u < 1$, the orthant probability for $\Sigma$ is $\frac{6S}{\pi^3}$ where $S=$

$$\begin{array}{ll} - \: 4 \: V( u, \frac{-a}{\sqrt{d}}, \frac{u}{\sqrt{d}}, \frac{-2 u}{\sqrt{3 b}}, -\sqrt{\frac{a}{3b}} ) & - \: 2 \: V( u, \frac{-a}{\sqrt{d}}, \frac{u^2}{\sqrt{3 d}}, \frac{2\sqrt{2}}{3}, -\sqrt{\frac{a}{3f}} ) \\ - \: 2 \: V( \frac{1}{2}, \frac{-u}{2}, \sqrt{a}, \frac{4 u^2}{\sqrt{i}}, \frac{a}{\sqrt{3 i}}) & - \: 2 \: V( \frac{1}{2}, \frac{-u}{2}, \frac{m}{2 \sqrt{h}}, \frac{l}{2}\sqrt{\frac{a}{h i}}, u^2 m \sqrt{\frac{2}{f i}} ) \\ - \: 4 \: V( \frac{1}{2}, \frac{-3u}{2 \sqrt{d}}, \frac{a}{\sqrt{d q}}, -2 u^2 \sqrt{\frac{3}{b q}}, -\sqrt{\frac{a}{3 b}} ) & - \: 2 \: V( \frac{1}{2}, \frac{-3u}{2 \sqrt{d}}, \frac{a}{2 \sqrt{d h}}, \sqrt{\frac{6}{h}}, -\sqrt{\frac{a}{3 f}} ) \\ - \: 2 \: V( \frac{u}{2}, \frac{-\sqrt{3}}{2}, \frac{-a}{2 \sqrt{b}}, \frac{l u}{2 \sqrt{b p}}, -u \sqrt{\frac{6 a}{f p}} ) & - \: 2 \: V( \frac{u}{2}, \frac{-j}2, u \sqrt{\frac{a}{n}}, \frac{4 j u}{\sqrt{i n}}, \frac{a}{\sqrt{3 i} }) \\ - \: 2 \: V( \frac{u}{2}, \frac{-j}2,\frac{u^2 m}{2 \sqrt{3 b}}, \frac{j l}{2}\sqrt{\frac{a}{3 b i}}, u^2 m \sqrt{\frac{2}{f i}} ) & + \: 2 \: V( u, \frac{-a}{\sqrt{d}}, \frac{u^2}{\sqrt{3 d}}, \frac{l}{3 \sqrt{b}}, u^2 \sqrt{\frac{2 a}{3 b f}} ) \\ + \: 2 \: V( u, \frac{-a}{\sqrt{d}}, u^2 \sqrt{\frac{3}{d}}, \frac{a}{\sqrt{b}}, u^2 \sqrt{\frac{2}{b}} ) & + \: 2 \: V( \frac{1}{2}, \frac{-u}{2}, \frac{\sqrt{a}}{2}, \frac{\sqrt{3}}{2}, -u^2 \sqrt{\frac{2}{g}} ) \\ + \: 2 \: V( \frac{1}{2}, \frac{-u}{2}, \frac{\sqrt{a}}{2}, 2 \sqrt{\frac{2}{k}}, \frac{-a}{\sqrt{3 g k}}) & + \: 2 \: V( \frac{1}{2}, \frac{-u}{2}, \frac{m}{2 \sqrt{h}}, \sqrt{\frac{2 a}{h k}}, \frac{-m}{\sqrt{k f} }) \\ + \: 2 \: V( \frac{1}{2}, \frac{-3u}{2 \sqrt{d}}, \frac{1}{2}\sqrt{\frac{a}{d}}, \frac{1}{2}\sqrt{\frac{3 a}{b}}, u^2 \sqrt{\frac{2}{b}} ) & + \: 2 \: V( \frac{1}{2}, \frac{-3u}{2 \sqrt{d}}, \frac{a}{2 \sqrt{d h}}, \frac{l}{2}\sqrt{\frac{3}{b h}}, u^2 \sqrt{\frac{2 a}{3 b f}} ) \\ + \: 4 \: V( \frac{1}{2}, \frac{-1}{2}, \frac{-2 u}{\sqrt{q}}, \frac{1}{2}\sqrt{\frac{a}{q}}, \frac{-1}{2} ) & + \: 4 \: V( \frac{u}{2}, \frac{-1}{2}, -\sqrt{\frac{2}{3}}, -\sqrt{\frac{a}{3 f}}, -u \sqrt{\frac{2}{f}} ) \\ + \: 4 \: V( \frac{u}{2}, \frac{-1}{2}, \frac{-3}{2}\sqrt{\frac{a}{e}}, -2 u^2 \sqrt{\frac{2}{e f}}, -a \sqrt{\frac{3}{f g}} ) & + \: 4 \: V( \frac{u}{2}, \frac{-1}{2}, \frac{-3}{2}\sqrt{\frac{a}{e}}, \frac{-u}{2 \sqrt{e}}, -u \sqrt{\frac{6}{g}} ) \\ + \: 2 \: V( \frac{u}{2}, \frac{-\sqrt{3}}{2}, \frac{-1}{2}\sqrt{\frac{a}{d}}, \frac{-3u}{2}\sqrt{\frac{a}{d p}}, u \sqrt{\frac{2}{p}} ) & + \: 2 \: V( \frac{u}{2}, \frac{-\sqrt{3}}{2}, \frac{-a}{2 \sqrt{b}}, u^2 \sqrt{\frac{2}{b}}, \sqrt{\frac{3 a}{f}} ) \\ + \: 2 \: V( \frac{u}{2}, -\sqrt{a}, \frac{-u^2}{\sqrt{e}}, -2 \sqrt{\frac{2 a}{3 e}}, \frac{-1}{3}\sqrt{\frac{a}{g}} ) & + \: 2 \: V( \frac{u}{2}, -\sqrt{a}, \frac{-u^2}{\sqrt{e}}, -a \sqrt{\frac{3}{e}}, -2 u^2 \sqrt{{2}{g}} ) \\ + \: 2 \: V( \frac{u}{2}, -\sqrt{a}, \frac{-u}{\sqrt{n}}, -8 u \sqrt{\frac{a}{3 g n}}, \frac{1}{3}\sqrt{\frac{a}{g}} ) & + \: 2 \: V( \frac{u}{2}, -\sqrt{a}, -u^2 \sqrt{\frac{3}{d}}, -a \sqrt{\frac{a}{d g}}, 2 u^2 \sqrt{\frac{2}{g}} ) \\ + \: 2 \: V( \frac{u}{2}, \frac{-j}{2}, \frac{u^2}{2}\sqrt{\frac{a}{e}}, \frac{j}{2}\sqrt{\frac{3}{e}}, -u^2 \sqrt{\frac{2}{g}} ) & + \: 2 \: V( \frac{u}{2}, \frac{-j}{2}, \frac{u^2}{2}\sqrt{\frac{a}{e}}, 2 j \sqrt{\frac{2}{e k}}, \frac{-a}{\sqrt{3 g k} }) \\ + \: 2 \: V( \frac{u}{2}, \frac{-j}{2}, \frac{u^2 m}{2 \sqrt{3 b}}, j \sqrt{\frac{2 a}{3 b k}}, \frac{-m}{\sqrt{k f} }) & . \end{array}$$

For instance, for $u=\frac{1}{2}$, we have $\frac{6S}{\pi^3} \approx 0.11502383599812541648615657162020006115$.