# What is the probability that a multivariate Gaussian random variable is greater than zero?

I am looking for a way to find the probability that $p(x > 0)$, where the vector $x$ has a multivariate Gaussian distribution $$x = \begin{bmatrix} x_1 \\ x_2 \end{bmatrix} \sim \mathcal{N}\left(\begin{bmatrix} \mu_1 \\ \mu_2 \end{bmatrix}, \begin{bmatrix} \Sigma_{11} & \Sigma_{12} \\ \Sigma_{21} & \Sigma_{22} \end{bmatrix}\right).$$ So I want to find $p(x_1 > 0, x_2 > 0)$. $x_1$ and $x_2$ are both scalars.

If $\Sigma_{12} = 0$, then this is easy, because then $x_1$ and $x_2$ are independent, meaning $$p(x_1 > 0, x_2 > 0) = p(x_1 > 0)p(x_2 > 0) = \left(\frac{1}{2}\mbox{erf}\left(\frac{\mu_1}{\sqrt{2\Sigma_{11}}}\right) + \frac{1}{2}\right)\left(\frac{1}{2}\mbox{erf}\left(\frac{\mu_2}{\sqrt{2\Sigma_{22}}}\right) + \frac{1}{2}\right).$$ However, I can't figure out how to do this for a non-independent Gaussian distribution. Of course approximating it numerically (just taking tons of samples and checking them) is an option, but I need something computationally more efficient. Any suggestions are appreciated.

Edit:

So I'll outline what I've done so far. First, for simplicity, I've assumed $\mu_1 = \mu_2 = 0$. (Adding $\mu$ later shouldn't be too difficult.) To simplify the problem further, I also assume that all parameters are scaled to give unit variance. So, $$x = \begin{bmatrix} x_1 \\ x_2 \end{bmatrix} \sim \mathcal{N}\left(\begin{bmatrix} 0 \\ 0 \end{bmatrix}, \begin{bmatrix} 1 & \rho \\ \rho & 1 \end{bmatrix}\right),$$ with $\rho$ the correlation coefficient. Then, I consider the problem of $p(x \leq 0)$ instead of $p(x \geq 0)$. This is exactly the same problem, but the notation will be slightly simpler. Next, for easy notation, I've defined $$\phi(x) = \int_{-\infty}^x \exp\left(-\frac{1}{2}s^2\right) \, ds = \sqrt{\frac{\pi}{2}} \left(1 + \mbox{erf}\left(\frac{x}{\sqrt{2}}\right)\right),$$ which will be useful to have. In particular, we now know that $$c\phi\left(\frac{x}{c}\right) = \int_{-\infty}^x \exp\left(-\frac{1}{2} \frac{s^2}{c^2}\right) \, ds,$$ $$\frac{d\phi\left(\frac{x}{c}\right)}{dx} = \frac{1}{c} \exp\left(-\frac{1}{2}\frac{x^2}{c^2}\right),$$ $$\int_{-\infty}^{x} \phi\left(\frac{s}{c}\right) \, ds = x\phi\left(\frac{x}{c}\right) + 2c\exp\left(-\frac{1}{2}\frac{x^2}{c^2}\right).$$ Now we can get to solving the problem. We can write the probability (in the general case) as $$p(x \leq 0) = \int_{-\infty}^0 \int_{-\infty}^0 \frac{1}{\sqrt{\left|2\pi\begin{bmatrix} \Sigma_{11} & \Sigma_{12} \\ \Sigma_{21} & \Sigma_{22} \end{bmatrix}\right|}} \exp\left(-\frac{1}{2} \begin{bmatrix} x_1 \\ x_2 \end{bmatrix}^T \begin{bmatrix} \Sigma_{11} & \Sigma_{12} \\ \Sigma_{21} & \Sigma_{22} \end{bmatrix}^{-1} \begin{bmatrix} x_1 \\ x_2 \end{bmatrix}\right) \, dx_1 \, dx_2.$$ If we write $\Delta = \Sigma_{11} - \Sigma_{12}\Sigma_{22}^{-1}\Sigma_{21}$, then this becomes (through the matrix inversion lemma and the matrix determinant lemma) $$p(x \leq 0) = \int_{-\infty}^0 \int_{-\infty}^0 \frac{1}{\sqrt{\left|2\pi\Delta\right|\left|2\pi\Sigma_{22}\right|}} \exp\left(-\frac{1}{2} (x_1 - \Sigma_{12}\Sigma_{22}^{-1}x_2)^T \Delta^{-1} (x_1 - \Sigma_{12}\Sigma_{22}^{-1}x_2) - \frac{1}{2} x_2^T \Sigma_{22}^{-1} x_2\right) \, dx_1 \, dx_2.$$ Now we will use the fact that all quantities are in fact scalars, to simplify notation. In addition, we will also use $\Sigma_{11} = \Sigma_{22} = 1$ and $\Sigma_{12} = \Sigma_{21} = \rho$. This gives us $$p(x \leq 0) = \int_{-\infty}^0 \int_{-\infty}^0 \frac{1}{2\pi\sqrt{1 - \rho^2}} \exp\left(-\frac{1}{2} \frac{(x_1 - \rho x_2)^2}{1 - \rho^2}\right) \exp\left(-\frac{1}{2} x_2^2\right) \, dx_1 \, dx_2.$$ Changing the integration limits of the inner integral gives us $$p(x \leq 0) = \int_{-\infty}^0 \int_{-\infty}^{-\rho x_2} \frac{1}{2\pi\sqrt{1 - \rho^2}} \exp\left(-\frac{1}{2} \frac{x_1^2}{1 - \rho^2}\right) \exp\left(-\frac{1}{2} x_2^2\right) \, dx_1 \, dx_2.$$ Solving the inner integral turns the above into $$p(x \leq 0) = \int_{-\infty}^0 \frac{1}{2\pi} \phi\left(\frac{-\rho x_2}{\sqrt{1 - \rho^2}}\right) \exp\left(-\frac{1}{2} x_2^2\right) \, dx_2.$$ This is the difficult integral to solve. We can apply integration by parts to it, which would result in $$p(x \leq 0) = \left[\frac{1}{2\pi} \phi\left(\frac{-\rho x_2}{\sqrt{1 - \rho^2}}\right) \phi\left(x_2\right)\right]_{-\infty}^0 + \int_{-\infty}^0 \frac{1}{2\pi} \frac{\rho}{\sqrt{1 - \rho^2}} \exp\left(-\frac{1}{2}\frac{\rho^2 x_2^2}{1 - \rho^2}\right) \phi\left(x_2\right) \, dx_2.$$ Substituting $x_2 = \frac{\sqrt{1 - \rho^2}}{\rho} y$ (assuming $\rho > 0$) will turn the above into $$p(x \leq 0) = \frac{\left(\phi(0)\right)^2}{2\pi} + \int_{-\infty}^0 \frac{1}{2\pi} \exp\left(-\frac{1}{2} y^2\right) \phi\left(\frac{\sqrt{1 - \rho^2}}{\rho} y\right) \, dy,$$ which is where we get stuck.

Effectively, we have turned the problem of solving the integral $$\int_{-\infty}^0 \frac{1}{2\pi} \phi\left(\frac{-\rho}{\sqrt{1 - \rho^2}} x\right) \exp\left(-\frac{1}{2} x^2\right) \, dx$$ into that of solving the integral $$\int_{-\infty}^0 \frac{1}{2\pi} \phi\left(\frac{\sqrt{1 - \rho^2}}{\rho} y\right) \exp\left(-\frac{1}{2} y^2\right) \, dy,$$ which is just as hard to solve.

Edit

Just now I found out, through some strange coincidence, that the solution for the simplified problem (with mean zero) is $$p(x \leq 0) = \frac{1}{4} + \frac{\sin^{-1}(\rho)}{2\pi}.$$ However, incorporating the mean $\mu$ into this relation seems to be impossible. So the problems remains yet unsolved.

Edit

I found out in literature that there is no closed-form solution. So I'll just stick with the Matlab mvncdf function.

the first step is to invert the covariance matrix so that you can express the joint distribution function as $$f(x_1,f_2) = \frac{1}{N} e^{-(A_{11} x_1^2 + A_{12}x_1 x_2+ A_{22}x_2^2)}$$ where $N$ is a normalization constant.

Then by integrating from $x_1 = 0$ to infinity you find the conditional distribution of $x_2$ given that $x_1 > 0$ will be product of an exponential of a quadratic in $x_2$ and an erf also involving $x_2$.

Now you want to integrate that for $x_2 = 0$ to infinity. The integral can be done by parts using the fact that $$\int \text{erf }(z) = z \text{ erf }(z) + \frac{e^{-z^2}}{\sqrt{\pi}}$$

• Thanks for the suggestions! They got me a step further. Sadly, integration by parts results in an integral of the same form, but with slightly different coefficients, so this does not result in a solution. I've added my work so far in the question statement. I believe it is what you suggested? – Hildo Bijl Apr 24 '15 at 13:33