Marginal distribution of Gaussians 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?
 A: There is an easy way to do this kind of problem and a hard way. 
The easy way is to notice that 
\begin{align}
y & = f + \varepsilon \\
 & \varepsilon \sim \mathcal{N}(0, \sigma^2I) 
\end{align} 
The expectation and variance of $ y $ give 
\begin{align}
E[y] & = E[f] + E[\varepsilon] = 0 \\
V[y] & = V[f] + V[\varepsilon]  = K + \sigma^2 I
\end{align}
Sums of normals are normal, so $ y \sim \mathcal{N}(0, K + \sigma^2 I) $
The hard way is to multiply out the pdf's and refactor it. Inside the exponential term, 
\begin{align}
y^T I/\sigma^2 y - 2 f^T  I/\sigma^2 y + f^T I/\sigma^2 f + f^T K^{-1} f \\
\end{align}
Move terms outside of the exponential that do not have to do with $f$ and complete the quadratic form, 
\begin{align}
& f^T A^{-1} f - 2f^Tb \\
A &= [ K^{-1} + I/\sigma^2]^{-1} \\
b & = I/\sigma^2 y \\
\end{align}
The completed quadratic form will have mean $ h $ defined as 
\begin{align}
h & = A b \\
c & =  - b^T A b \\
(x & - h) A^{-1} (x-h) + c \\
\end{align} 
Pull out the $ c $ term and notice that this is Gaussian, so it has a known normalizing constant. The $ c $ combines with the terms involving $ y $ that we moved outside of the integral in the first step, and find another Gaussian:
\begin{align}
y ^T[ I/\sigma^2 + I/\sigma^2[I/\sigma^2 + K^{-1} ]^{-1} I/\sigma^2] y
\end{align}
The Woodbury matrix inversion formula gives another form of the covariance matrix, 
$$
I/\sigma^2 + I/\sigma^2[I/\sigma^2 + K^{-1} ]^{-1} I/\sigma^2 = \sigma^2I + K 
$$
See https://en.wikipedia.org/wiki/Woodbury_matrix_identity and let 
\begin{align}
A & = \sigma^2I \\
U & = I \\
V & = I \\
C & = K^{-1}
\end{align}
Therefore,
$$ y \sim \mathcal{N} (0, \sigma^2I + K )$$
