in a post from over four years ago, Przemo gave the following formula for the integral over a Gaussian function over the positive reals in three dimensions (denoted here as $\mathbb{R}^3_+$) with inverse covariance $A=\left(\begin{array}{rrr}a & a_{12} & a_{13}\\a_{12}& b&a_{23}\\a_{13}&a_{23}&c\end{array}\right)$:
\begin{multline} \int_{\mathbb{R}^3_+} \mathrm{d}x \; \exp\left(-\frac{1}{2}x^\top A x\right) = \\ \frac{\pi}{\sqrt{2 \det A}}\int\limits_0^\infty \mathrm{erfc}\left(u \sqrt{\frac{-2 a_{12}a_{13}a_{23}+a_{13}^2 b+a_{12}^2 c}{\det(A)}}\right)e^{-u^2}\cdot \\ \left[-\mathrm{erfc}\left(\sqrt{a} \frac{a_{12} c-a_{13}a_{23}}{a_{13} \sqrt{\det(A)}} u\right)+\mathrm{erfc}\left(\sqrt{a} \frac{a_{12}a_{23}-a_{13} b}{a_{12} \sqrt{\det(A)}} u\right)\right] du \end{multline}
This integral can be done using the techniques outlined here, yielding a closed-form solution to the original problem.
I would like to understand the derivation of the intermediate result given above. In his original post, Przemo details the calculation in two dimensions: complete the square in the quadratic form; do the first integral, which yields an error function; expand the error function and do the remaining integral term by term. How can you adapt this approach to 3D?
Steps undertaken so far:
Completing the square in one variable, say $x$, leaves me with $$\int_{\mathbb{R}_+^2} \mathrm{d}y\mathrm{d}z \exp\left(-\frac{1}{2} \frac{\mathrm{det}\,A_3}{\mathrm{det}\,A_2}z^2\right) \exp\left(-\frac{1}{2} \frac{\mathrm{det}\, A_2}{a}(y-m z)^2\right) \left[1 - \mathrm{erf}\left(\frac{a_{12}y+a_{13}z}{\sqrt{2a}}\right) \right] $$ where $A_2=\begin{pmatrix} a & a_{12}\\ & b\end{pmatrix}$, $A_3$ as defined above, and $m$ is a function of the coefficients of the matrices. However, I do not know how to proceed from there: expanding the error function to do the integral in y, say, is a nightmare due to the constant term in z; I also did not find a way to do a coordinate transform à la $s=a_{12}y+a_{13}z$ or something similar.
Indeed, the formula above looks more like one completed the square in two of the variables independently; but what happened to the cross-term? I cannot find a factorisation of the exponent that would allow me to complete two integrals over the half-line with only one variable left in the error function yielded by the integral.
Any help / hint would be greatly appreciated! Thank you in advance for your time :)
