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.

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.

1 Answer
1

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

$\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 BijlApr 24 '15 at 13:33