# Normal orthant probability in six variables

Let $$\mathbf{X} \sim N(\mathbf{0}, \mathbf{\Sigma})$$ be a $$6$$-dimensional Gaussian vector with covariance matrix of the form $$\mathbf{\Sigma} = \begin{pmatrix} 1 & c \\ c & 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 & c & c/2 & c/2 \\ 1/2 & 1 & 1/2 & c/2 & c & c/2 \\ 1/2 & 1/2 & 1 & c/2 & c/2 & c \\ c & c/2 & c/2 & 1 & 1/2 & 1/2 \\ c/2 & c & c/2 & 1/2 & 1 & 1/2 \\ c/2 & c/2 & c & 1/2 & 1/2 & 1 \end{pmatrix}$$ Here $$c$$ is some constant between $$0$$ and $$1$$.

I am interested in $$\Pr(\mathbf{X} \geq \mathbf{0})$$, i.e. the probability that all of the six coordinates are non-negative/positive. 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, 3, 4$$ dimensions. The analysis seems very tedious and nasty though, and I am not sure if the same techniques generalize easily to $$6$$ dimensions.

For my example of $$6$$ dimensions, my search has not yet returned any literature, attempting to solve this problem even in special cases with additional structure like above. I also tried computing this probability with Mathematica, but it cannot solve it analytically, and even numerically it seems to have a hard time to return exact results.

My question is: is there any way to find the orthant probability for such a structured $$6$$-dimensional matrix? Is there any literature I am missing? Or can someone solve this analytically?

Edit: The paper https://ieeexplore.ieee.org/document/1054159 describes a derivation of the orthant probability for four variables, where \begin{align} \mathbf{\Sigma} = \begin{pmatrix} 1 & c \\ c & 1 \end{pmatrix} \otimes \begin{pmatrix} 1 & d \\ d & 1 \end{pmatrix}. \end{align} Perhaps it is possible to use similar techniques for this case as well (using a similar path of integration approach), although the reduction at the bottom of page 389 would not go from a $$4$$-dimensional orthant probability to a $$2$$-dimensional one, which can readily be evaluated, but from $$6$$ dimensions to $$4$$ dimensions, and adding another integral over these orthant probabilities.

Still, this reference might be useful -- it might inspire similar techniques/approaches to this problem.

• In case this helps: this can alternatively be phrased as the probability that two 3D vectors $X, Y$, each with individual correlation matrix $\begin{pmatrix} 1 & 1/2 & 1/2 \\ 1/2 & 1 & 1/2 \\ 1/2 & 1/2 & 1 \end{pmatrix}$ between their three coordinates, and with correlation $c$ between the two vectors, both lie in the positive orthant.
– TMM
Commented Jun 7, 2019 at 12:09
• I was thinking along the lines of cholesky decomposition of $\Sigma$ to get the matrix $L$ such that $y=Lx$ is ~ $N(0,I)$. Conditions of $x_i \geq 0$ will map to a different set of linear inequations involving $y_i$s. At this point the probabilities could be calculated via simulation as $y_i$s are independent standard normals. But I couldn't think of a better analytical approach. Commented Jun 12, 2019 at 2:06
• @DebNandy Using an R library for numerically evaluating orthant probabilities, I can already get reasonably accurate numerical results, but I was wondering if an analytical approach exists for solving this problem. (Analytically, that decomposition just shifts the problem from the correlated random variables to dependent boundaries for the integrals.)
– TMM
Commented Jun 12, 2019 at 18:40
• I posted my (attempt at an) answer to a similar question here math.stackexchange.com/a/4301412/16878. It won't give nice-looking formulas (unless there are some simplifications) but it should (hopefully) allow one to calculate the result exactly for any $n$. Unless, of course, there are errors and it actually turn out to be wrong. :) Commented Nov 9, 2021 at 16:54

This is not a complete answer, but maybe this will help you to get a formula in terms of the distribution function of the multivariate normal distribution.

Lets denote by $$X = (X_1,\ldots,X_6) \sim N(0,\Sigma)$$ the original random vector. Set $$I = \{1,2,3,4,5,6\}$$. For $$J \subset I$$ I will denote by $$F_J$$ the distribution function of the random variables $$X_{J_1},\ldots,X_{J_{m}}$$, where $$m = |J|$$. For example taking $$J = \{1,3,5\}$$ we have $$m = 3$$ and $$F_J(x) = \Pr[X_1 \leq x_1, X_3 \leq x_2, X_5 \leq x_3]$$ denotes the distribution function of $$(X_1,X_3,X_5)$$ at $$x=(x_1,x_2,x_3)$$ which is again jointly normal. The joint distribution function of $$X$$ is then $$F_I$$ and for $$J=\{2,6\}$$ we have $$F_J(x) = \Pr[X_2 \leq x_1, X_6 \leq x_2]$$ and so forth.

Your orthant probability can be thought of as computing the measure of a box under the joint probability measure and is linked to the inclusion / exclusion principle. To express this in terms of the joint distribution function I will follow Chapter 10 in Nelsen (2006): An introduction to copulas.

For two vectors $$a < b$$ (meaning $$a_i < b_i$$ for all $$i$$) the corresponding box is simply $$B = \times_{i=1}^6 [a_i,b_i]$$, and $$c=(c_1,\ldots,c_6)$$ will denote a vertex (or corner) of $$B$$, that is: $$c_k$$ is either equal to $$a_k$$ or $$b_k$$. By equation (2.10.1) in Nelsen (2006) the measure of $$B$$ under $$V_{F_I}$$ (the measure induced by $$F_I$$) is now given as \begin{align} V_{F_I}(B) = \sum \text{sgn}(c)F_I(c), \end{align} where the sum is taken over all vertices $$c$$ of $$B$$, and $$\text{sgn}(c)$$ is given as \begin{align} \text{sgn}(c) = 1, \text{if } c_k = a_k \text{ for an even number of k's and -1 otherwise.} \end{align} Here the sum will have $$2^6 = 64$$ terms which seems unmanageable, but in your case two simplifications are possible.

First, your orthant is a box with $$a = (0,\ldots,0)$$ and $$b=(\infty,\ldots,\infty)$$. If a corner $$c$$ has an infinity as a component, the joint distribution $$F_I(c)$$ will actually reduce to a lower dimensional one. For example in case of $$c=(0,\infty,0,\infty,0,\infty)$$ we have $$F_I(c) = F_J(0,0,0)$$ where $$J=\{1,3,5\}$$.

Second, with the specific covariance structure that you have in mind these lower dimensional distributions might be the same for some index sets, meaning that it is possible that $$F_{J_1} = F_{J_2}$$ for some $$J_1$$ and $$J_2$$ with $$|J_1|=|J_2|$$.

What you would need to check now is which lower dimensional distributions are the same (this depends on $$\Sigma$$) and which sign they get. Then you can collect equal terms in the sum and rewrite the sums in terms of the dimension of the lower dimensional distributions going from one to six.

In your case you will have for each dimension multiple terms with different coefficients, but if your random vector had the same covariance between all components you would end up with six terms only (and the coefficients would be related to the binomial coefficients).

I hope this helps or gives maybe a different direction to tackle your problem. Good luck!

• By Chapter 10 I suppose you mean Section 2.10? Thanks for the reference - I never even heard of copulas before. The inclusion-exclusion sounds very similar to a paper of David, academic.oup.com/biomet/article-abstract/40/3-4/458/…, which however explains that for even dimensions, this reduction does not work; we get $V_H(B) = H(\infty, \dots, \infty) \pm \dots + H(0, \dots, 0)$ with a plus-sign for $H(0, \dots, 0) = V_H(B)$ by symmetry. Subtracting $V_H(B)$ from both sides, this only gives us an identity in the other quantities, but does not solve $V_H(B)$.
– TMM
Commented Jun 13, 2019 at 23:24