Marginalizing multivariate Gaussian distribution 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.
 A: 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 :)
