2
$\begingroup$

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

$\endgroup$
3
$\begingroup$

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

$\endgroup$
  • 4
    $\begingroup$ An even simpler way is to replace Probability with NProbability (adds one additional character to the OP's code). $\endgroup$ – DumpsterDoofus Jan 24 '15 at 0:41
  • $\begingroup$ Yep... Simpler. $\endgroup$ – David G. Stork Jan 24 '15 at 0:42
  • $\begingroup$ Nice, thank you very much $\endgroup$ – Filanto Renika Jan 24 '15 at 0:52
  • $\begingroup$ The general integral for the OPs $n=3$ case does actually have a closed-form solution: posted below. $\endgroup$ – wolfies Jan 24 '15 at 16:47
9
$\begingroup$

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.

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

Your Answer

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

Not the answer you're looking for? Browse other questions tagged or ask your own question.