While working through the exercises in [Mathematics for machine learning][1] 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][2] and [some of the answers in this site][3].

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][4]:
$$
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.


  [1]: https://mml-book.github.io/
  [2]: https://en.wikipedia.org/wiki/Multivariate_normal_distribution#Marginal_distributions
  [3]: https://stats.stackexchange.com/a/401398/110833
  [4]: http://www.maths.qmul.ac.uk/~ig/MTH5118/Notes11-09.pdf