3

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

J. M.'s missing motivation
  • 124,525
  • 11
  • 401
  • 574
Filanto Renika
  • 33
  • 3

2 Answers2

9

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.

wolfies
  • 8,722
  • 1
  • 25
  • 54
  • 1
    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. –  Jan 24 '15 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. – wolfies Jan 24 '15 at 16:59
  • Not quite all done. Let's help the OP by giving him the closed form solution for his particular covariance matrix. – David G. Stork Jan 24 '15 at 19:21
  • 2
    @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$. – Michael E2 Jan 24 '15 at 23:22
3

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

Mr.Wizard
  • 271,378
  • 34
  • 587
  • 1,371
David G. Stork
  • 41,180
  • 3
  • 34
  • 96