1
$\begingroup$

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.

$\endgroup$
1
$\begingroup$

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}} $$

$\endgroup$
  • $\begingroup$ 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? $\endgroup$ – Hildo Bijl Apr 24 '15 at 13:33

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.