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 $\matthcal{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.

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.