2
$\begingroup$

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

$\endgroup$
  • $\begingroup$ @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?) $\endgroup$ – fedja Jun 17 '19 at 21:55
  • $\begingroup$ @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? $\endgroup$ – fedja Jun 17 '19 at 22:08
  • $\begingroup$ @TMM In other words, the original expression was $1/4$ above the truth. OK, now it makes sense. :) $\endgroup$ – fedja Jun 18 '19 at 0:36
  • $\begingroup$ 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. $\endgroup$ – TMM Jun 18 '19 at 0:38

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service, privacy policy and cookie policy

Browse other questions tagged or ask your own question.