Marginalizing multivariate Gaussian distribution

Question

While working through the exercises in Mathematics for machine learning I have encountered a claim (Eq. (6.68)) that the marginal of a two-dimensional normal distribution $\mathcal{N}(x, y |\mathbf{\mu}, \mathbf{\Sigma})$ is simply $\mathcal{N}(x |\mu_x, \Sigma_{xx})$. This claim is echoed by Wikipedia and some of the answers in this site.

Yet, the direction calculation gives a different result: $$ p(x,y) = \frac{1}{2\pi\sqrt{\Sigma_{xx}\Sigma_{yy} - \Sigma_{xy}^2}} e^{-\frac{1}{2}\left[\Sigma_{xx}(x-\mu_x)^2 + \Sigma_{yy}(y-\mu_y)^2 + 2\Sigma_{xy}(x-\mu_x)(y-\mu_y)\right]},\\ p(x) = \int dy p(x,y) = \frac{1}{\sqrt{2\pi}\left(\Sigma_{xx} - \frac{\Sigma_{xy}^2}{\Sigma_{yy}}\right)} e^{-\frac{1}{2}\left(\Sigma_{xx} - \frac{\Sigma_{xy}^2}{\Sigma_{yy}}\right)(x-\mu_x)^2}, $$ where I use $$ \int dy e^{-py^2 + qy} = \sqrt{\frac{\pi}{p}}e^{\frac{q^2}{4p}}, $$ which is just a shortcut for completing the square.

My result seems to be confirmed by a different representation of the bivariate Gaussian distribution: $$ p(x,y) = \frac{1}{2\pi\sigma_1\sigma_2\sqrt{1-\rho^2}} e^{\frac{-1}{2(1-\rho^2)}\left[\left(\frac{x-\mu_1}{\sigma_1}\right)^2 -2\rho \left(\frac{x-\mu_1}{\sigma_1}\right) \left(\frac{y-\mu_2}{\sigma_2}\right) + \left(\frac{y-\mu_2}{\sigma_2}\right)^2\right]}, $$ where, as notation implies, $$ \Sigma_{xx} - \frac{\Sigma_{xy}^2}{\Sigma_{yy}} = \frac{1}{\sigma_1^2}. $$

It is hard to believe in a widespread error of this magnitude. On the other hand, the calculation is very straightforward... I will appreciate clarifications.

Answer 1

I think that a few results from Robb J. Muirehead, Aspects of Multivariate Statistic Theory, might be useful.

Definition 1.2.3. The $m\times 1$ random vector $\mathbf{X}$ is said to have an $m$-variate normal distribution if, for every $\mathbf{a}\in\mathbb{R}^m$, the distribution of $\mathbf{a}^T\mathbf{X}$ is univariate normal.

Theorem 1.2.6. If $\mathbf{X}$ is $N_m(\boldsymbol{\mu},\boldsymbol{\Sigma})$ and $\mathbf{B}$ is $k\times m$, $\mathbf{b}$ is $k\times 1$, then $$\mathbf{Y}=\mathbf{BX}+\mathbf{b}\quad\text{is}\quad N_k(\mathbf{B}\boldsymbol{\mu}+\mathbf{b}, \mathbf{B}\boldsymbol{\Sigma}\mathbf{B}^T)$$ Proof. Direct consequence of Definition 1.2.3, and of a few well known facts.

Theorem 1.2.7. If $\mathbf{X}$ is $N_m(\boldsymbol{\mu},\boldsymbol{\Sigma})$ then the marginal distribution of any subset of $k$ ($<m$) components of $\mathbf{X}$ is $k$-variate normal.
Proof. This follows directly from the Definition, of from Theorem 1.2.6. For example $$\mathbf{X}=\begin{bmatrix} \mathbf{X}_1 \\ \mathbf{X}_2 \end{bmatrix},\quad \boldsymbol{\mu}=\begin{bmatrix}\boldsymbol{\mu}_1\\\boldsymbol{\mu}_2\end{bmatrix},\quad \boldsymbol{\Sigma}=\begin{bmatrix}\boldsymbol{\Sigma}_{11} & \boldsymbol{\Sigma}_{12}\\ \boldsymbol{\Sigma}_{21} & \boldsymbol{\Sigma}_{22}\end{bmatrix}$$ Putting $$\underset{k\times m}{\mathbf{B}}=[\mathbf{I}_k:0],\qquad\mathbf{b}=\mathbf{0}$$ in Theorem 1.2.6 shows immediately that $\mathbf{X}_1\sim N(\boldsymbol{\mu}_1,\boldsymbol{\Sigma}_{11})$.

This is much clearer and easier, IMHO. But I love linear algebra :)