Find closed form for quadruple integral I am trying to find a closed form of the following integral
$$
\int _0^{\infty }\int _0^x\int _0^y\int _0^z \exp \left( -\frac{a x^2}{2}-\frac{b y^2}{2}-\frac{c z^2}{2}-\frac{d w^2}{2} \right) \,\mathrm{d}w\,\mathrm{d}z\,\mathrm{d}y\,\mathrm{d}x
$$
where $a,b,c,d>0$ are some constants.
My idea is to change variable to to polar system by letting
$$
x=\frac{r \cos (\alpha ) \cos (\beta ) \cos (\theta )}{\sqrt{a}}, \quad
y=\frac{r \cos (\alpha ) \cos (\beta ) \sin (\theta )}{\sqrt{a}}
$$
$$
z=\frac{r \sin (\alpha ) \cos (\beta )}{\sqrt{c}}, \quad w=\frac{r \sin (\beta )}{\sqrt{d}}
$$
This reduces the original integral into
$$
\int_0^{\tan ^{-1}\left(\sqrt{\frac{b}{a}}\right)} \frac{\sin (\theta ) \tan ^{-1}\left(\sin (\theta ) \sqrt{\frac{d}{b+c \sin ^2(\theta )}}\right)}{\sqrt{a b d \left(b+c \sin ^2(\theta )\right)}} \, d\theta
$$
But then I get stuck here.

PS: I am interested in this because I found that
$$
\int _0^{\infty }\int _0^x \exp \left(-\frac{a x^2}{2}-\frac{b y^2}{2}\right)\,\mathrm{d}y\,\mathrm{d}x
=
\frac{\tan ^{-1}\left(\sqrt{\frac{b}{a}}\right)}{\sqrt{a b}}
$$
and
$$
\int _0^{\infty }\int _0^x\int _0^y \exp \left( -\frac{a x^2}{2}-\frac{b y^2}{2}-\frac{c z^2}{2} \right) \,\mathrm{d}z\,\mathrm{d}y\,\mathrm{d}x
=
\frac{\sqrt{\pi/2 } }{\sqrt{a b c}}
\left(\tan ^{-1}\left(\sqrt{\frac{c}{b}}\right)-\tan ^{-1}\left(\sqrt{\frac{a c}{b (a+b+c)}}\right)\right)
$$
So I am trying to generalize this. Maybe this is already known?
 A: This integral is up to normalization constant integral of multivariate gaussian distribution. Due too lack of cross terms we get that there are 4 independent zero mean gaussian variables involved.
We can begin by rewriting this a probability of an event.
Let $X,Y,Z,W$ independent normally distributed with zero mean and possibly different variances.
Then integral reduces to:
$$ \mathbb{P} (0 < W < Z < Y < X) $$
Which in turn is the same as:
$$ \mathbb{P}(W > 0 \land Z-W >0 \land Y-Z >0 \land X-Y >0) $$
and one can observe that $W,Z-W,Y-Z,X-Y$ are correlated jointly normal variables.
Probability that all components of a jointly normal vector are positive is called orthant probability and in general doesn't have closed form expression. I suppose this is quite special case with covariance matrix almost diagonal so maybe there are some articles about how to approach this special case.
For the case with 3 or less variables formulas are known (cf this question for example) and I guess they would coincide with what you've found.
