Proof of the Derivation of the marginal and conditional Gaussian Given a marginal Gaussian distribution for x and a conditional Gaussian distribution for y given x in the form 
$$p(x) = N(x|\mu, \Lambda^{-1})$$
$$p(y|x) = N(y|Ax + b, L^{-1})$$
the marginal distribution of y and the conditional distribution of x given y are given by 
$$p(y) = N(y|A\mu +b, L^{-1} + A\Lambda^{-1}A^T)$$
$$p(x|y) = N(x|\Sigma{A^TL(y-b) + \Lambda \mu}, \Sigma)$$
where $\Sigma = (\Lambda + A^TLA)^{-1}$.
This is taken from the textbook Pattern Recognition and Machine Learning.
The textbook has a proof but it is very brief and I could not fully understand the proof. Does anyone have a simple and detailed proof that can help me understand this fact? Thanks.
 A: Let's write the RV $X$, $Y$ as
$$
X = \mu + \varepsilon _x \\
Y = AX + b + \varepsilon _y
$$
with $ \varepsilon _x \sim \mathcal N (0, \Lambda ^{-1})$, $ \varepsilon _y \sim \mathcal N (0, L ^{-1})$. Now pluging in X in the second equation above gives
$$
Y = A\mu + b + A\varepsilon _x + \varepsilon _y.
$$
This is a linear combination of normal distributed random variables and as such itself normal distributed with expectation $A\mu +b$ and covariance matrix $A\Lambda ^{-1}A^T + L^{-1}$ ($var(AX)=Avar(X)A^T$). From this you get $p(y)$.
The second fact is a bit more complicated and involves some tedious calculation. From Bayes Theorem it follows that
$$
p(x|y) \propto p(y|x)p(x) \\
\propto \exp ((y-Ax-b)^TL(y-Ax-b) + (x-\mu)^T\Lambda(x-\mu)).
$$
If you multiply everything out and then factor out $x$ you get to
something proportional to
$$
\exp( (x-(\Lambda + A^TLA)^{-1}A^TL(y-b))^T(\Lambda + A^TLA)(x-(\Lambda + A^TLA)^{-1}AL(y-b)))
$$
which is proportional to your given normal.
A: I would like to also give a algebraic response to this question showing that by multiplying the Gaussians together and integrating you can get this result. I did this using the covariance matrix defined as $L, \Lambda$ instead of $L^{-1}$ and $\Lambda^{-1}$ because that is my preference.  
\begin{align}
y|x \sim &  N(Ax + b, L) \\
x \sim & N(\mu, \Lambda) 
\end{align}
\begin{multline}
y^T L^{-1} y - 2y^T L^{-1}(Ax + b) + (Ax+b)^T L^{-1} (Ax + b) + \\ x^T\Lambda^{-1} x - 
2x^T\Lambda^{-1} \mu + \mu^T \Lambda^{-1} \mu = \\
y^T L^{-1} y- 2y^T L^{-1} Ax - 2y^T L^{-1} b  + x^T A^T L^{-1} A x + 2x^T A^T L^{-1} 
b + b^T L^{-1} b \\  x^T\Lambda^{-1} x - 2x^T\Lambda^{-1} \mu + \mu^T \Lambda^{-1} \mu
\end{multline}  
As before make quadratic forms out of the $ x $ terms to integrate them out:
\begin{align}
x^T ( \Lambda^{-1} + A^T L^{-1} A) x - 2x^T ( \Lambda^{-1} \mu + A^T L^{-1} y - A^T L^{-1} b) + 
\mu^T \Lambda^{-1} \mu + b^T L^{-1} b
\end{align}
\begin{align}
V &= [A^TL^{-1} A + \Lambda^{-1} ]^{-1} \\
u &= [ \Lambda^{-1} \mu + A^T L^{-1} y - A^T L^{-1} b]\\
h & = Vu \\
h^TV^{-1}h &= u^T V u \\
c & = \mu^T \Lambda^{-1} \mu + b^T L^{-1} b \\
\end{align}
$ x $ is Gaussian and integrates to its normalizing constant. The term $ c- K 
$, $ K = u^TVu $ factors into the $ y $ exponential giving, 
\begin{align}
y^T L^{-1} y -2y^TL^{-1} b + c - K 
\end{align}
And, 
\begin{multline}
u^T V u = \\
 [ \Lambda^{-1} \mu + A^T L^{-1} y - A^T L^{-1} b]^T V  [ \Lambda^{-1} \mu + A^T L^{-1} y - A^T L^{-1} b] = \\
\mu^T \Lambda^{-1} V \Lambda^{-1} \mu + y^T L^{-1} A V A^T L^{-1} y + b^T L^{-1} AV A^T L^{-1} b 
\\ 
+ 2y^T L^{-1} A V \Lambda^{-1} \mu - 2\mu^T \Lambda^{-1} V A^T L^{-1} b - 2y^T L^{-1} A V 
A^T L^{-1} b 
\end{multline}
Together this becomes:
\begin{multline}
y^T [ L^{-1} - L^{-1} A V A^T L^{-1} ] y - 2y^T [L^{-1} b + L^{-1} A V \Lambda^{-1} \mu - L^{-1} A V A^T L^{-1} b] + \\ \mu^T[ \Lambda^{-1} - \Lambda^{-1} V 
\Lambda^{-1} ]\mu + b^T[L^{-1} - L^{-1} 
A V 
A^T L^{-1}] b - 2 b^T L^{-1} A V \Lambda^{-1} \mu 
\end{multline}
Which simplifies slightly, 
\begin{equation}
y^T R^{-1} y - 2y^T g + b^T R^{-1} b +  \mu^T[ \Lambda^{-1} - \Lambda^{-1} V\Lambda^{-1}]\mu - 2 b^T 
L^{-1} A V \Lambda^{-1} \mu
\end{equation}
\begin{align}
R^{-1} = L^{-1} - L^{-1} A V A^T L^{-1}  = [L + A \Lambda A^T]^{-1}   \\
g = [L^{-1} b + L^{-1} A V \Lambda^{-1} \mu - L^{-1} A V A^T L^{-1} b]   
\end{align}
Furthermore, $ Rg  = A\mu +b $:
\begin{align}
R g &= R [L^{-1} b + L^{-1} A V \Lambda^{-1} \mu - L^{-1} A V A^T L^{-1} b] \\
&=  R[ (L^{-1} - L^{-1} A VA^T L^{-1}) b + L^{-1} A V \Lambda^{-1} \mu  ] \label{amu1} \\
& = b + R L^{-1} A V \Lambda^{-1} \mu \label{amu2}
\end{align}
After canceling out the $ b $ term on both sides I get I solve for the other 
term, equating coefficients:
\begin{align}
R L^{-1} A V \Lambda \mu = & [L + A \Lambda A^T] L^{-1} A V \Lambda^{-1} \mu = A \mu 
\label{amu3} \\
= & A  + A\Lambda A^T L^{-1} A = A [\Lambda^{-1} + A^T L^{-1} A]  \Lambda \label{amu4} \\
& \text{(using the definition of $ V $)} \notag \\
\end{align}
\begin{align}
I + \Lambda A^T L^{-1} A = & I + A^T L^{-1} A \Lambda \label{amu5}\\
\Lambda A^T L^{-1} A = & \Lambda A^T L^{-1} A  \label{amu6}
\end{align}
And therefore $ Rg = A\mu + b $ and $ R^{-1} A \mu = L^{-1} AV \Lambda^{-1} \mu. $
Simplifying equation the equation for $y$ I get, 
\begin{equation}
y^T R^{-1} y - 2y^TR^{-1} (A\mu + b)+ b^TR^{-1} b  + \mu^T[ \Lambda^{-1} - \Lambda^{-1} V 
\Lambda^{-1} ]\mu - 2 b^T L^{-1} A V \Lambda^{-1} \mu.
\end{equation}
Notice,
$$
\Lambda - \Lambda^{-1} V \Lambda^{-1}   = A^T R^{-1} A
$$
To show this:
\begin{align}
A^T R^{-1} A & =  \Lambda^{-1} - \Lambda^{-1} V \Lambda^{-1} \\
& = \Lambda^{-1} -  \Lambda^{-1} I ( I \Lambda I + F^{-1}) I \Lambda^{-1}  
\end{align}
$ F^{-1} = A^TL^{-1} A $,
\begin{align}
\Lambda^{-1} - \Lambda^{-1} I ( I \Lambda^{-1} I + F^{-1}) I \Lambda^{-1}  & = [\Lambda + F] ^{-1} 
\end{align}
By the Woodbury identity. 
\begin{align}
A^T R^{-1} A & =  [\Lambda + F] ^{-1}\\
A^{-1} R A^{-1} & = \Lambda + F \\
R & = A \Lambda A^T + L 
\end{align}
Since both sides are equal, the `sandwich' term is equal to $ A^T R^{-1} A $
 I get the 
following:
\begin{equation}
y^T R^{-1} y - 2y^TR^{-1} (A\mu + b)+ b^TR^{-1} b  + \mu^T A^T R^{-1} A \mu - 2 
b^T R^{-1} 
A\mu.
\end{equation}
Which gives the results that $ y $ is distributed normally. 
$$
y \sim \mathcal{N} (A\mu + b, R)
$$
The other part is much easier, and the algebra is not very tedious. The first the long equation given for $ y $ can be simplified when conditioning on $x$ since everything not involving $x$ can be subsumed into the constant of proportionality. This leaves only, 
\begin{equation}
-2x^T A^T L^{-1} y + x^T A^T L^{-1} x + 2 x^T A^T L^{-1} b + X^T \Lambda^{-1} x - 2x^T \Lambda ^{-1} \mu 
\end{equation}
Then refactoring this term, 
\begin{equation}
x^T [ A^T L^{-1} A + \Lambda^{-1} ] x - 2x^T [A^T L^{-1} (y - b) + \Lambda^{-1} \mu] + K 
\end{equation} Where $ K $ is some constant.
Call $[ A^T L^{-1} A + \Lambda^{-1} ] = \Sigma^{-1}$. Therefore, 
\begin{equation}
x \sim \mathcal{N} ( \Sigma [ A^T L^{-1}(y - b) + \Lambda^{-1} \mu], \Sigma)
\end{equation}
Remember I defined my covariance matrix to the inverse of what was given in the question.
