Probability of multivariate normal being positive on each coordinate

Question

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).

Answer 1

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, ∞}]

yield the answer

(* 1/3 *)

Answer 2

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.