I want to calculate the log marginal likelihood for a Gaussian Process regression, for that and by GP definition I have the prior: $$ p(\textbf{f} \mid X) = \mathcal{N}(\textbf{0} , K)$$
Where $ K $ is the covariance matrix given by the kernel. And the likelihood is (a factorized gaussian):
$$ p(\textbf{y} \mid \textbf{f}, X) = \mathcal{N}(\textbf{f} ,\sigma_n²I)$$
where $ \textbf{f} $ are the training outputs (the values of the function) with some random gaussian noise term with variance $ \sigma_n²I $ So the log marginal likelihood is calculated as follows: $$ p(\textbf{y}\mid X) = \int p(\textbf{y} \mid \textbf{f}, X) p(\textbf{f} \mid X) d\textbf{f}$$
which ends up in the solution
$$ y \sim \mathcal{N}(\textbf{0} , K + \sigma_n²I) $$ so the log marginal likehood is given by: $$ \log p(\textbf{y} \mid X) = -\frac{1}{2} \textbf{y}^T (K + \sigma_n²I)^{-1}\textbf{y} -\frac{1}{2} \log |K+\sigma_n²I| -\frac{n}{2}\log(2\pi)$$ ($n$ is the number of training points).
The book I'm currently following (Rasmussen, C. E., & Williams, C. K. (2006). Gaussian processes for machine learning, Page 19 ) uses Gaussian identities to get to the solution, specifically the the multiplication of two gaussian functions, but I'm a bit loss on how to apply them to get to the solution, ¿how do I calculate the marginal distribution?