# Probability of multivariate normal being positive on each coordinate

How can I find the probability that each coordinate of a specified multivariate normal distribution is positive? I tried the following, which I believed should work

mu = {0, 0, 0};
sigma = {{2, 1, 1}, {1, 2, 1}, {1, 1, 2}};
Probability[
x > 0 && y > 0 && z > 0, {x, y, z} \[Distributed]
MultinormalDistribution[mu, sigma]]


Unfortunately, for the output I just get the last line from the input (with mu and sigma replaced by their actual values). I don't see where the problem could possibly be since the matrix is positive definite. If I replace it by the identity matrix everything works fine (i.e. the output is 1/8).

The general integral does not have a closed form solution, so use NIntegrate:

mu = {0, 0, 0};
sigma = {{2, 1, 1}, {1, 2, 1}, {1, 1, 2}};
NIntegrate[
PDF[MultinormalDistribution[mu, sigma], {x, y, z}],
{x, 0, ∞}, {y, 0, ∞}, {z, 0, ∞}]


(* 0.25 *)

Check:

mu = {0, 0, 0};
sigma2 = {{1, 0, 0}, {0, 1, 0}, {0, 0, 1}};
NIntegrate[
PDF[MultinormalDistribution[mu, sigma2], {x, y, z}],
{x, 0, ∞}, {y, 0, ∞}, {z, 0, ∞}]


(* 0.125 *)

Note that symbolic integration and symbolic probability work for the two-dimensional version of this problem:

mu = {0, 0};
sigma = {{2, 1}, {1, 2}};
Probability[x > 0 && y > 0, {x, y} \[Distributed] MultinormalDistribution[mu, sigma]]


and

Integrate[PDF[MultinormalDistribution[mu, sigma], {x, y}], {x, 0, ∞}, {y, 0, ∞}]


(* 1/3 *)

• An even simpler way is to replace Probability with NProbability (adds one additional character to the OP's code). Commented Jan 24, 2015 at 0:41
• Yep... Simpler. Commented Jan 24, 2015 at 0:42
• Nice, thank you very much Commented Jan 24, 2015 at 0:52
• The general integral for the OPs $n=3$ case does actually have a closed-form solution: posted below. Commented Jan 24, 2015 at 16:47

The probability that the OP seeks is known as the multivariate Normal orthant probability. Correctly, for the $n=3$ cased posed here, the general integral DOES in fact have a closed -form solution, though Mma cannot (currently) obtain it.

In particular, given a zero mean vector and variance-covariance matrix:

$$\Sigma =\left( \begin{array}{ccc} 1 & \rho _{\text{xy}} & \rho _{\text{xz}} \\ \rho _{\text{xy}} & 1 & \rho _{\text{yz}} \\ \rho _{\text{xz}} & \rho _{\text{yz}} & 1 \\ \end{array} \right)$$

$\dots$ the standardised trivariate Normal orthant probability is:

$$P(X>0,Y>0,Z>0) \quad = \quad \frac{1}{8} + \frac{\text{ArcSin}\left[\rho _{\text{xy}}\right]+ \text{ArcSin}\left[\rho _{\text{xz}}\right]+ \text{ArcSin}\left[\rho _{\text{yz}}\right]}{4 \pi }$$

For application and more detail, see, for instance:

• Chapter 6 of our book: Rose and Smith, Mathematical Statistics with Mathematica (Section 6.4C) $\rightarrow$ a free download is available at: http://www.mathstatica.com/book/bookcontents.html , or

• Stuart and Ord (1994), Kendall's Advanced Theory of Statistics (6th edition): section 15.10 and 15.11.

Example

Let $(X,Y,Z)$ have a standardised multivariate Normal with zero mean vector, and variance covariance matrix:

sigma = {{1, 27/34, 22/23}, {27/34, 1, 4/5}, {22/23, 4/5, 1}}


Then, the exact orthant probability is given immediately as:

P3 = 1/8 + (ArcSin[27/34] + ArcSin[22/23] + ArcSin[4/5])/(4 Pi)


... which, to 10 decimal places, is:

N[P3, 10]
0.3732564868


By contrast, the approach using numerical integration can be unreliable here:

NIntegrate[ PDF[MultinormalDistribution[{0, 0, 0}, sigma], {x, y, z}], {x, 0, Infinity}, {y, 0, Infinity}, {z, 0, Infinity}]


NIntegrate::slwcon: Numerical integration converging too slowly ...

0.371907

Numerical integration can sometimes be awkward, unreliable, or slow (as in this example) ... and having an exact closed-form instantaneous solution is a better way to proceed, if possible.

From general to standardised

Given a zero mean vector, what if our variance-covariance matrix is not in a standardised form (with 1's along the main diagonal)? We can easily convert it into standardised form. If, say, our variance-covariance matrix is:

  sigma = {{3, 1/3, 3/2}, {1/3, 2/3, -1}, {3/2, -1, 4}}


... then the standardised variance-covariance matrix S is:

 A = DiagonalMatrix[Diagonal[sigma]^(-(1/2))];
S = A.sigma.A


All done.

• I eagerly looked up your book because I wanted to see "more detail"... but there isn't any more detail, just a restatement of the same formula and a reference to Stuart and Ord's book.
– user484
Commented Jan 24, 2015 at 16:51
• The book also gives the closed-form solution for the bivariate case, and then uses these exact formulae as benchmarks to check the performance/accuracy of Mma's numerical integration routines. In case you are about to eagerly consult Stuart and Ord, I should perhaps caution that you won't find very much more detail there either for the $n=2$ or 3 cases as to results (though more as to derivation)... but they do also provide useful references for higher dimensions. Commented Jan 24, 2015 at 16:59
• Not quite all done. Let's help the OP by giving him the closed form solution for his particular covariance matrix. Commented Jan 24, 2015 at 19:21
• @Rahul I believe the result can be obtained by a linear change of variables that diagonalizes and normalizes the quadratic form of the pdf (to Exp[k (u^2+v^2+w^2)]). The orthant is transformed to a domain that is subtended by a spherical triangle with angles $\pi/2+ \sin^{-1}\rho$, where $\rho$ runs over the nondiagonal coefficients of $\Sigma$. The result is the proportion of the sphere that is cut out (by the symmetry of the integral). A similar result is true for $n=2$. Commented Jan 24, 2015 at 23:22