In this post series we study ellipsoidal cores of Multivariate Normal Distributions [MVNs]. We defined a “core” as the volume enclosed by a selected contour surface of constant probability density. We constructed cut-off distributions by setting the probability density to zero outside the core. In previous posts we have already discussed volume integrals which would give us a relation between the covariance matrix of the infinitely extended MVN and the covariance matrix of the core. We will evaluate all required integrals and factors for the 3-dimensional case in this post.
In this post we evaluate all of the required integrals by applying an “integration by parts” and by using the so called error-function. The integrals we get confronted with in 3 dimensions, are more complex than their counterparts in two dimensions. We will later see that this pattern continues – our integration problem is in general worse for an odd number of dimensions than for an even number.
A 2-dimensional example showed us that we must take a normalization factor for the cut-off distribution into account. In this post we evaluate a respective integral over a MVN’s probability density up to a certain Mahalanobis distance analytically. In a forthcoming post we will instead handle the integration for the normalization factor by an evaluation of the cumulative chi-square function in n dimensions.
While we cover the theory in this post, we will look at experimental data for the 3D case in the next post. The results will give us confidence in a generalization pattern which will allow us to cover higher numbers of dimensions (n ≥ 4).
Previous posts
Prerequisistes: We use abbreviations and conventions defined in the 1st post of this series. Note that some authors use the abbreviations BVD/MVD instead of our acronyms BVN/ MVN, respectively. I assume that you are familiar with the creation of a MVN for a random vector YN (with correlations between its components) by applying an invertible, linear transformation to a special MVN for a standardized random vector Z, composed of independent Gaussians. Z has a simple probability density with spherical contour surfaces. I.e., a geeral MVN YN can be generated from Z by applying some matrix MY:
\[ \pmb{Y}_N \,= \, \pmb{\operatorname{M}}_Y \, \pmb{Z} \,. \]
MY can also be understood as a matrix which mediates a coordinate transformation. See post 2 for details.
Integrals to solve for the n-dimensional case
From post 2 we know that we have to solve the following volume integral to get the covariance matrix CovC of the contents of a MVN’s ellipsoidal core.
\[
\pmb{\operatorname{Cov_C}} \,=\, \pmb{\operatorname{M}}_Y \left[ {1\over (2\,\pi)^{n/2} } \, \underset{ V_{n,D}^{sph}}{\iint …\int} \, \pmb{Z} \pmb{Z}^T \, \operatorname{exp} \left( -\, {1\over 2} \, \pmb{Z}^T \, \pmb{Z} \,\right) \, dz_1\,dz_2 \, ….\, dz_n \right] \pmb{\operatorname{M}}_Y^T \,. \tag{1}
\]
The term in the square brackets actually defines a matrix, because Z=(Z1, Z2, …, Zn)T. The integration (2) has to be performed over the volume of a n-dimensional sphere of radius r=D :
\[ \begin{align}
{\small{\text{Volume of n-dim sphere of radius }}} d_m=D \,:\, & \quad V_{n,D}^{sph} \,, \\[10pt]
{\small{\text{Volume of n-dim sphere of radius }}} d_m=r \,:\, & \quad V_{n,r}^{sph} \,, \\[10pt]
{\small{\text{Surface of n-dim sphere of radius }}} d_m=r \,:\, & \quad A_{n,r}^{sph} \,.
\end{align} \]
Note: Regarding dimensionality subscripts, I deviate a bit from other authors, who prefer to use (n-1) as a subscript of surfaces to indicate the lower dimensionality of the manifold in the ℝn.
See post 2 for reasons why we can work with a spherical volume although the core of YN is encapsulated in an ellipsoidal surface. The spherical volume is just the result of mapping the ellipsoidal core of YN ‘s MVN by a linear (coordinate) transformation [MY]-1 to Z and its spherical probability density. MY fulfills the following relations regarding YN ‘s covariance matrix CovY = ΣY:
\[ \begin{align} \pmb{\operatorname{\Sigma}}_{Y} \,&=\, \pmb{\operatorname{Cov_Y}} \,= \, \pmb{\operatorname{M}}_{\small Y}\, \pmb{\operatorname{M}}_{\small Y}^T \,, \tag{2} \\[10pt]
\pmb{\operatorname{\Sigma}}_{Y}^{-1} \,&=\, \left(\pmb{\operatorname{M}}_{\small Y}^T\right)^{-1} \pmb{\operatorname{M}}_{\small Y}^{-1} \,=\, \left(\pmb{\operatorname{M}}_{\small Y}^{-1}\right)^T \pmb{\operatorname{M}}_{\small Y}^{-1} \,. \tag{3}
\end{align} \, .
\]
CovY = ΣY is a symmetric, invertible and positive definite matrix.
Three dimensions, only …
For three dimensions, the exponential in integral (2) gives us is the probability density of the standardized 3-dim MVN of a standardized random vector Z3d
\[ \pmb{Z}_{3d} \,= \, \begin{pmatrix} Z_1 \\ Z_2 \\ Z_3 \end{pmatrix} \,, \tag{4} \]
with no correlation between its three independent components
\[ \begin{align}
& \text{standardized 3-dim MVN without correlation of the marginals: } \\[10pt]
& \pmb{Z}_{3d} \, \sim \pmb{\mathcal{N}}_3 \left( \pmb{0}, \, \pmb{\operatorname{I}} \right), \,\,\, {\small with} \,\,\, Z_1, Z_2, Z_3 \,\sim\, \mathcal{N}_1 \left( 0, \,1 \right) \,. \tag{5}
\end{align} \]
Z3d has a simple probability density function g3,z(z), characterizing the density at the end-point of a vector z3d = (z1, z2, z3)T:
\[ \begin{align}
& g_{3, {\small Z}}(\pmb{z}_{3d}) \,=\, {1 \over (2\pi) } \, e^{ – \, {\Large 1 \over \Large 2} \, { \Large |\pmb{z}_{3d}|}^2 } \,. \tag{6} \\[10pt]
& \pmb{z}_{3d} \, = \, (z_1, \, z_2, \, z_3)^T\,, \,\,\, |\pmb{z}_{3d}| \,=\, r \,=\,\sqrt{z_1^2 + z_2^2 + z_3^2 \, } \,. \tag{7}
\end{align} \]
A 3-dim target MVN Y3d can be created with the help of a (3×3)-matrix M3
\[\pmb{Y}_{3d} \,=\, \pmb{\operatorname{M}}_3 \, \pmb{Z}_{3d} \,. \tag{8} \]
We name Y3d ‘s cov-matrix Σ3 = Cov3,Y
\[ \pmb{\Sigma}_3 \,= \, \pmb{\operatorname{Cov}}_ {3, \pmb{\operatorname{Y}}} \,=\, \pmb{\operatorname{M}}_3\, \pmb{\operatorname{M}}_3^T \,. \tag{9}
\]
Our integral becomes
\[
\pmb{\operatorname{Cov}}_{3,\pmb{\operatorname{C}}} \, =\, \pmb{\operatorname{M}}_3 \left[ {1\over (2\,\pi)^{3/2} } \, \underset{ V_{3,D}^{sph}}{\iint } \, \pmb{Z}_{3d} \pmb{Z}_{3d}^T \, \operatorname{exp} \left( -\, {1\over 2} \, \pmb{Z}_{3d}^T \, \pmb{Z}_{3d} \,\right) \, dz_1\,dz_2 \,dz_3 \right] \pmb{\operatorname{M}}_3^T \,. \tag{10}
\]
The “volume” V3,Dsph of the spherical core (for dm=D) in (10) actually is the volume of a 3-sphere. To make the notation look more familiar, we replace z1, z2, z3 by x, y and z:
\[ \pmb{z}_{3d} \, = \, (x, \, y, \, z )^T\,, \,\,\,\,|\pmb{z}_{3d}| \,=\, r \,=\,\sqrt{x^2 + y^2 + z^2\, } \,. \tag{11} \]
Our integral then can be written as
\[
\pmb{\operatorname{Cov}}_{3,\pmb{\operatorname{C}}} \, =\, \pmb{\operatorname{M}}_3 \left[ {1\over (2\,\pi)^{3/2} } \, \underset{ V_{3,D}^{sph}}{\iint …\int} \, (x,\, y,\, z)^T\circ (x, \,y,\, z) \, \, \operatorname{exp} \left( -\, {1\over 2} \, r^2 \right) \, dx \,dy \,dz \, \right] \pmb{\operatorname{M}}_3^T \,. \tag{12}
\]
Due to the matrix product of a line vector with a column vector we get matrix elements, which have to investigate more closely.
Integrals for the off-diagonal elements
Let us write down a first example for the off-diagonal elements of the matrix Cov3,C. As we already saw in post 2, they all should become zero. Let us look e.g. at the (x,y) combination:
\[ \operatorname{Int}_{xy,D} \,=\, \left[ {1\over (2\,\pi)^{3/2} } \, \underset{ V_{3,D}^{sph}}{\iint} \, x \, y * \operatorname{exp} \left( -\, {1\over 2} \, r^2 \,\right) \, dx \,dy \,dz \right] \,. \tag{13} \]
Switching to polar coordinates (r, θ, φ)
\[ \begin{align} x \,&=\, r * \sin \theta \, \cos( \varphi) \,, \\ y \,&=\, r * \sin \theta \, \sin(\varphi) \,, \\ z \,&=\, r * \cos \theta \,. \end{align} \tag{14} \]
we indeed find
\[ \begin{align} \operatorname{Int}_{xy,D} \,&=\, {1\over (2\,\pi)^{3/2} } \, {\int_0^D \int_0^{\pi} \int_0^{2\pi} } \, r^2 \, \sin^2 \theta \, \cos\varphi \, \sin\varphi \,\, \operatorname{exp} \left( -\, {1\over 2} \, r^2 \,\right)\,\, r^2\,\sin\theta \, d \varphi \, d\theta \, dr \\[10pt]
&=\, {1\over 2\,\pi)^{3/2} } \, \int_0^D \, r^4 \, \operatorname{exp} \left( -\, {1\over 2} \, r^2 \right) \, dr \, \left[ \int_0^{\pi} \sin^3 \theta \, d\theta \, \right] \, \int_0^{2\pi} \, \, \cos\varphi \, \sin\varphi \, d\varphi \\[10pt]
& =\, 0\,.
\end{align} \tag{15} \]
Or
\[ \begin{align} \operatorname{Int}_{yz,D} \,&=\, {1\over 2\,\pi)^{3/2} } \, {\int_0^D \int_0^{\pi} \int_0^{2\pi} } \, r^2 \, \sin \theta \, \cos \theta \, \cos\varphi \, \sin\varphi \,\, \operatorname{exp} \left( -\, {1\over 2} \, r^2 \,\right)\,\, r^2\,\sin\theta \, d \varphi \, d\theta \, dr \\[10pt]
&=\, {1\over 2\,\pi)^{3/2} } \, \int_0^D \, r^4 \, \operatorname{exp} \left( -\, {1\over 2} \, r^2 \, \right) \, dr \, \left[ \int_0^{\pi} \sin^2 \theta \cos\theta \, \, d\theta \, \right] \, \int_0^{2\pi} \, \sin\varphi \, d\varphi \\[10pt]
& =\, 0\,.
\end{align} \tag{16} \]
Analogously for the (x,z)-related term. I.e.: CovC in our 3-dim case becomes a diagonal matrix W3,diag :
\[ \pmb{\operatorname{Cov}}_ {3, \pmb{\operatorname{C}}} \,=\, \pmb{\operatorname{W}}_{3,diag} \,=\, \begin{pmatrix}
W_{xx} & 0 & 0 \\ 0 & W_{yy} & 0 \\ 0 & 0 & W_{zz} \end{pmatrix} \, . \tag{17} \]
Based on elementary symmetry arguments, we had already predicted this in post 2 – even for the general n-dimensional case.
Diagonal matrix elements
We can already conclude by symmetry reasons that Wxx = Wyy = Wzz. But let us show this explicitly:
\[ \begin{align} W_{xx} \,&=\, {1\over 2\,\pi)^{3/2} } \, {\int_0^D \int_0^{\pi} \int_0^{2\pi} } \, r^2 \, \sin^2 \theta \, \cos^2(\varphi) \,\, \operatorname{exp} \left( -\, {1\over 2} \, r^2 \,\right)\, r^2 \, \sin\theta \, d \varphi \, d\theta \, dr \\[10pt]
&=\, {1\over 2\,\pi)^{3/2} } \, {1\over 2}\, \int_0^{D} r^4 \, \operatorname{exp} \left( -\, {1\over 2} \, r^2 \,\right) \, dr \, \left[ \int_0^{\pi} \, \sin^3 \theta \, d\theta \, \right] \, \int_0^{2\pi} \cos^2(\varphi) \, d \varphi \, \\[10pt]
&=\, {4 \, \pi \over 3} *{1\over (2\,\pi)^{3/2} } \, {1\over 2}\, \int_0^{D} r^4 \, \operatorname{exp} \left( -\, {1\over 2} \, r^2 \,\right) \, \, dr \, \,.
\end{align} \tag{18} \]
\[ \begin{align} W_{yy} \,&=\, {1\over (2\,\pi)^{3/2} } \, {\int_0^D \int_0^{\pi} \int_0^{2\pi} } \, r^2 \, \sin^2 \theta \, \sin^2(\varphi) \,\, \operatorname{exp} \left( -\, {1\over 2} \, r^2 \,\right)\, r^2 \, \sin\theta \, d \varphi \, d\theta \, dr \\[10pt]
&=\, {1\over (2\,\pi)^{3/2} } \, {1\over 2}\, \int_0^{D} r^4 \, \operatorname{exp} \left( -\, {1\over 2} \, r^2 \,\right) \, \, dr \, \left[ \int_0^{\pi} \, \sin^3 \theta \, d\theta \, \right] \, \int_0^{2\pi} \sin^2(\varphi) \, d \varphi \, \\[10pt]
&=\, {4 \, \pi \over 3} *{1\over (2\,\pi)^{3/2} } \, {1\over 2}\, \int_0^{D} r^4 \, \operatorname{exp} \left( -\, {1\over 2} \, r^2 \,\right) \, \, dr \, \,.
\end{align} \tag{19} \]
\[ \begin{align} W_{zz} \,&=\, {1\over (2\,\pi)^{3/2} } \, {\int_0^D \int_0^{\pi} \int_0^{2\pi} } \, r^2 \, \cos^2 \theta \, \,\, \operatorname{exp} \left( -\, {1\over 2} \, r^2 \,\right)\, r^2 \, \sin\theta \, d \varphi \, d\theta \, dr \\[10pt]
&=\, {1\over (2\,\pi)^{3/2} } \, {1\over 2}\, \int_0^{D} r^4 \, \operatorname{exp} \left( -\, {1\over 2} \, r^2 \,\right) \, dr \, \left[ \int_0^{\pi} \, \cos^2 \theta \, \sin\theta \, d\theta \, \right] \, \int_0^{2\pi} \, d \varphi \, \\[10pt]
&=\, {4 \, \pi \over 3} *{1\over (2\,\pi)^{3/2} } \, {1\over 2}\, \int_0^{D} r^4 \, \operatorname{exp} \left( -\, {1\over 2} \, r^2 \,\right) \, \, dr \, \,.
\end{align} \tag{20} \]
So, indeed,
\[ W_{xx} \,=\, W_{yy} \,=\, W_{zz} \,=: \, W_3 \,. \tag{21} \]
This does not surprise us, as it just confirms conclusions based on elementary symmetry considerations in post 2. Its nevertheless nice to see that the different integrals just provide the right factors.
The radial part
We have to evaluate the remaining radial integral. It can be done with the help of an integration by parts. I go through the steps, because we will need their pattern for a generalization to higher dimensions in a forthcoming post. We make use of the following derivatives:
\[ \begin{align}
{\partial \over \partial x} \left[ x^3 \, e^{-1/2 \, x^2} \right] \,&=\, 3\, x^2 \, e^{-1/2 \, x^2} \,-\, \left( x^4\, e^{-1/2 \, x^2} \right) \,, \tag{22} \\[10pt]
{\partial \over \partial x} \left[ 3\, x \, e^{-1/2 \, x^2} \right] \,&=\, 3\, e^{-1/2 \, x^2} \,-\, \left(-\, 3 x^2\, e^{-1/2 \, x^2} \right) \,, \tag{23} \\[10pt]
{\partial \over \partial x} \left[ {\sqrt{\pi} \over \sqrt{2}} \, \operatorname{erf} \left( {x\over \sqrt{2}} \right) \right] \,&=\, e^{-1/2 \, x^2} \,. \tag{24} \\[10pt]
\end{align} \]
With the help of these relations we once again integrate by parts
\[ \begin{align}
\operatorname{I}_{rad} \, &= \, \int_0^{D} \, r^4 \, \operatorname{exp} \left( -\, {1\over 2} \, r^2 \,\right) \, \, dr \\[10pt]
& = \, \left[ – \, r^3 \, \operatorname{exp} \left( -\, {1\over 2} \, r^2 \right) \, \right]_0^D \,+\, \int_0^{D} \,3\,r^2 \, \operatorname{exp} \left( -\, {1\over 2} \, r^2 \,\right) \, \, dr \\[10pt]
&= \, \left[ – \,r^3 \, \operatorname{exp} \left( -\, {1\over 2} \, r^2 \right) \, \right]_0^D \,+ \, \left[ -\, 3\, r \, \operatorname{exp} \left( -\, {1\over 2} \, r^2 \right) \right]_0^D \,+ \, 3 \, \int_0^{D} \, \operatorname{exp} \left( -\, {1\over 2} \, r^2 \right) \, dr \\[10pt]
&=\, \left[ \, 3\, { \sqrt{\pi} \over \sqrt{2} }\, \operatorname{erf} \left(r/\sqrt{2}\right) \, – \, ( r^3 + 3\, r) \, \operatorname{exp} \left( -\, {1\over 2} \, r^2 \right) \, \right]_0^D \, .
\end{align} \tag{25} \]
You can find the antiderivative also in a list of integrals published by N.E. Korotkov and A.N. Korotkov; see [1]. “erf” is the error function
\[ \operatorname{erf}(x) \,=\, 2\, \Phi_0(x/\sqrt{2})\,; \quad \text{with} \,\,\, \Phi_0(y) \, = \, {1 \over \sqrt{2 \,\pi\, }} \, \int_0^y \, \exp (- {1 \over 2} \, t^2 ) \, dt \,. \tag{26} \]
Bringing the results together, we conclude
\[ \begin{align}
W_3 \, &= \, {4 \, \pi \over 3} *{1\over (2\,\pi)^{3/2} } \, {1\over 2}\, \operatorname{I}_{rad} \\[10pt]
&=\, {1 \over 3} \, {1 \over \sqrt{2\, \pi \,} } \, \left[ \, 3\, { \sqrt{\pi} \over \sqrt{2} }\, \operatorname{erf} \left(r/\sqrt{2}\right) \, – \, ( r^3 + 3\, r) \, \operatorname{exp} \left( -\, {1\over 2} \, r^2 \right) \, \right]_0^D \\[10pt]
&=\, {1 \over 2} \, \operatorname{erf} \left(D/\sqrt{2}\right) \,-\, {1 \over 3} \, {1 \over \sqrt{2\, \pi \,} } \, \left[ (D^3 + 3 \, D) \exp(- \,{1 \over 2} D^2) \right] \, .
\end{align} \tag{27} \]
Who said that life should be easy …
Normalization factor
When working with data samples we can numerically determine the covariance matrix Cov3,C for the data points enclosed in a given 3-dimensional core. From the 2-dimensional example (covered in the last post) we know already that we need a normalization factor to be able to estimate the full MVN’s covariance matrix Σ3 based on data for Cov3,C. This normalization factor is (by a mapping to the z-coordinates) given by
\[
F_{3,D} \,=\, {1 \over (2\pi)^{3/2} } \,*\, \underset{ V_{3,D}^{sph}} {\iiint} \, \operatorname{exp} \left( -\, {1\over 2} \, \pmb{z}^T \pmb{z} \,\right) \, dz_1\,dz_2 \,dz_3 \,. \tag{28}
\]
Switching to spherical polar coordinates we get:
\[ \begin{align}
F_{3,D} \,&=\, {1 \over (2\pi)^{3/2} } \,*\, \int_0^D \int_0^{\pi} \int_0^{2\pi} \, \operatorname{exp} \left( -\, {1\over 2} \, r^2 \,\right) \, r^2 \, \sin \theta \, dr \, d\theta \, d\varphi \\[10pt]
&=\, 2 * 2\,\pi \, {1 \over (2\pi)^{3/2} } \,*\, \int_0^D \, r^2 \operatorname{exp} \left( -\, {1\over 2} \, r^2 \,\right) \, dr \\[10pt]
&=\, 2 \, {1 \over \sqrt{2 \, \pi \,} } \,*\, \int_0^D \, r^2 \operatorname{exp} \left( -\, {1\over 2} \, r^2 \,\right) \, dr \,.
\end{align} \tag{29} \]
This is an already familiar integral. Picking up derivatives above, we integrate by parts again:
\[ \begin{align}
F_{3,D} \,&=\, {2 \over \sqrt{2 \, \pi \,} } \,*\ \left[ \, {\sqrt{\pi} \over \sqrt{2}} \, \operatorname{erf} ( {r \over \sqrt{2}} ) – r\, \operatorname{exp} \left( -\, {1\over 2} \, r^2 \,\right) \, \right]_0^D \\[10pt]
&=\, \operatorname{erf} ( {D \over \sqrt{2}} ) \,-\, {\sqrt{2} \over \sqrt{\pi}} \, D\, \operatorname{exp} \left( -\, {1\over 2} \, D^2 \,\right) \,.
\end{align} \tag{30} \]
The interested reader can verify that F3,D is identical to the cumulative chi-square function for 3 degrees of freedom. The probability density gχ2 of the chi-square distribution for 3 degrees of freedom is given by:
\[ g_{\chi^2}(x) \,=\, { x^{(n-2)/2} \, e^{-x/2} \over 2^{n/2} \, \Gamma(n/2) } \tag {31} \]
with Γ being the Gamma-function. To show the equivalence with (28) you must of course replace x by r2.
The relation between the core’s covariance matrix Cov3,C and the MVN’s Cov3,Y
From post 2 and post 3 we know that the relation between the core’s covariance matrix and the cov-matrix of the infinite MVN is given by
\[ \pmb{\operatorname{\Sigma}}_3 \,=\, \pmb{\operatorname{Cov}}_{3, \pmb{Y}} \,=\, { F_{3,D} \over W_3} \, \pmb{\operatorname{Cov}}_{3, \pmb{C}} \,. \tag{32} \]
I.e., we have succeeded with the evaluation all necessary terms and can verify the relation by numerical experiments.
Conclusion
Previous posts in this series have well prepared us to evaluate all integrals which determine the covariance matrix of a core of a MVN in 3 dimensions. The attentive reader has certainly detected patterns in the steps of a integration by parts, which he have used. These patterns will allow us later on to generalize the results for multidimensional MVNs in the ℝn. This in turn gives us a tool to analyze data samples in statistics and Machine Learning which resemble a MVN only inside a central dense region.
In the next post we will confirm the theoretical results of this post for a numerical test case in three dimensions.
Links and literature
[1] N.E. Korotkov and A.N. Korotkov, 2019, “Table of integrals related to error function” by Nikolai E. Korotkov, retired Leading Researcher, Voronezh Institute of Communications, Russia, edited by Alexander N. Korotkov, Professor, University of California, Riverside, USA
https://intra.engr.ucr.edu/~korotkov/papers/Korotkov-book-integrals.pdf