How to perform joint inference on multivariate normal variables? Suppose I have the following model:
$$\begin{aligned}
\text C &\sim \mathcal N \left(\mu, \delta^2\right) \\
\forall i: \text L_i | \text C = c &\sim \mathcal N \left(c, \lambda_i^2 \right) \\
\forall i: \text P_i | \text L_i = l_i &\sim \mathcal N \left(l_i, \sigma_i^2 \right)
\end{aligned}$$
In the above, $\text C$ and the $\text L_i$ are the hidden variables I'm reasoning about, and the $\text P_i$ are the observed variables that serve as evidence. Then:
$$\begin{aligned}
p \left(\text C, \vec {\text L} \middle | \vec {\text P} = \vec p \right) &\propto
p \left(\vec {\text P} = \vec p \middle | \text C, \vec {\text L} \right) p \left(\vec {\text L} \middle | \text C \right) 
p \left(\text C \right) \\
&=
p \left(\text C \right) \prod\limits_{i=1}^n 
p\left(\text P_i = p_i \middle| \text L_i\right) p\left(\text L_i \middle | \text C\right) \\
&=
\mathcal N \left(C \middle |\mu, \delta^2\right) \prod\limits_{i=1}^n \mathcal N \left(p_i \middle |\text L_i, \sigma_i^2 \right)
\mathcal N \left(\text L_i \middle| \text C, \lambda_i^2 \right) \\
&\propto
\mathcal N \left(C \middle |\mu, \delta^2\right)
 \prod\limits_{i=1}^n 
\mathcal N \left(\text L_i \middle| \frac {\lambda_i^{-2} \text C + \sigma_i^{-2} p_i}{\lambda_i^{-2} + \sigma_i^{-2}}, \left(\lambda_i^{-2} + \sigma_i^{-2} \right)^{-1} \right) \\
\end{aligned}$$
Presumably the above is a multivariate normal distribution, with some posterior mean $\vec \mu$ and some posterior covariance matrix $\mathbf \Sigma$. What are they?
A possible way to figure this out would be to describe the original problem vectorially, but I have no idea how to do that either. I think, since in the multivariate normal definition the mean vector seems to be defined as the vector of the means of each variable unconditionally on all others, we'd have:
$$\begin{aligned}
\text C, \vec {\text L} &\sim
\mathcal N \left(
\begin{bmatrix} \mu \\ \mu \\ \vdots \\ \mu\end{bmatrix},
\begin{bmatrix}
\delta^2 & \delta^2 & \cdots & \delta^2 \\
\delta^2 & \delta^2 + \lambda_1^2 & & \delta^2 \\
\vdots &&\ddots & \\
\delta^2 & \delta^2 & & \delta^2 + \lambda_n^2
\end{bmatrix}
\right)
\end{aligned}$$
Although I'm not sure the first row and column have the right values.
If the above is correct, then:
$$\begin{aligned}
\vec {\text P} | \vec {\text L} = \vec l
&\sim \mathcal N \left(\vec l,
\begin{bmatrix}
\sigma_1^2 && 0\\
& \ddots &\\
0 && \sigma_n^2
\end{bmatrix} \right) \\
\vec {\text L} | \text C = c
&\sim \mathcal N \left(
\begin{bmatrix}c \\ \vdots \\ c \end{bmatrix},
\begin{bmatrix}
\lambda_1^2 && 0 \\
& \ddots & \\
0 && \lambda_n^2
\end{bmatrix}
\right)
\end{aligned}$$
So for the inference step:
$$\begin{aligned}
p \left(\text C, \vec {\text L} \middle| \vec {\text P} = \vec p \right)
&\propto
\mathcal N \left(\text C \middle |\mu, \delta^2\right)
 \mathcal N \left(\vec p \middle | \vec {\text L},
\begin{bmatrix}
\sigma_1^2 && 0\\
& \ddots &\\
0 && \sigma_n^2
\end{bmatrix} \right) \mathcal N \left(\vec {\text L} \middle|
\begin{bmatrix}\text C \\ \vdots \\ \text C \end{bmatrix},
\begin{bmatrix}
\lambda_1^2 && 0 \\
& \ddots & \\
0 && \lambda_n^2
\end{bmatrix}
\right) \\
&=
\mathcal N \left(\text C \middle |\mu, \delta^2\right)
 \mathcal N \left(\vec p \middle | \vec {\text L}, \mathbf \Sigma \right) \mathcal N \left(\vec {\text L} \middle |\vec {\text C}, \mathbf \Lambda
\right) \\
&\propto
\mathcal N \left(\text C \middle |\mu, \delta^2\right)
 \mathcal N \left(\vec {\text L} \middle | \left( \mathbf \Sigma ^ {-1} + \mathbf \Lambda ^ {-1}\right) ^ {-1} \left( \mathbf \Sigma ^ {-1} \vec p + \mathbf \Lambda ^ {-1} \vec {\text C} \right), \left( \mathbf \Sigma ^ {-1} + \mathbf \Lambda ^ {-1} \right) ^ {-1} \right)
\end{aligned}$$
Except the above is just... the same thing, but written vectorially. So I'm stumped. What's next, how do I find the posterior mean and covariance matrix for the whole thing?
 A: Okay, so let me try to continue the vectorial thoughts. If I define $\vec {\text {CL}}$ as the joint vector of my latent variables, then there exists a matrix $\mathbf X$ such that $\vec L = \mathbf X \vec{\text {CL}}$:
$$\mathbf X = 
\begin{bmatrix}
0 & 1 & & 0 \\
\vdots &  & \ddots & \\
0 & 0 & & 1
\end{bmatrix}$$
Therefore:
$$\begin{aligned}
p \left(\vec {\text {CL}} \middle| \vec {\text P} = \vec p \right)
&= \frac {p \left(\vec {\text P} = \vec p \middle| \vec {\text {CL}} \right)p \left(\vec {\text {CL}} \right)}{p \left(\vec {\text P} = \vec p \right)}
\\ &\propto
\mathcal N \left(\vec p \middle | \mathbf X \vec{\text {CL}}, \mathbf \Sigma\right)
\mathcal N \left(\vec{\text {CL}} \middle |
\begin{bmatrix} \mu \\ \mu \\ \vdots \\ \mu\end{bmatrix},
\begin{bmatrix}
\delta^2 & \delta^2 & \cdots & \delta^2 \\
\delta^2 & \delta^2 + \lambda_1^2 & & \delta^2 \\
\vdots &&\ddots & \\
\delta^2 & \delta^2 & & \delta^2 + \lambda_n^2
\end{bmatrix}
\right)
\\ &=
\mathcal N \left(\vec p \middle | \mathbf X \vec{\text {CL}}, \mathbf \Sigma\right)
\mathcal N \left(\vec{\text {CL}} \middle | \vec \mu, \mathbf \Delta \right)
\\ &\propto
\exp {\left(\left( \vec p  - \mathbf X \vec{\text {CL}} \right)^\top \mathbf \Sigma^{-1} \left( \vec p  - \mathbf X \vec{\text {CL}} \right) +
\left(\vec{\text {CL}} - \vec\mu \right)^\top \mathbf \Delta^{-1} \left(\vec{\text {CL}} - \vec\mu \right)
 \right)}^{-\frac 1 2}
\\ &=
\exp {\left(
\begin{aligned}
&\left(\vec p ^\top - \vec{\text {CL}}^\top \mathbf X ^\top \right)  \left(\mathbf \Sigma^{-1} \vec p  - \mathbf \Sigma^{-1} \mathbf X \vec{\text {CL}} \right) + \\
&\left(\vec{\text {CL}}^\top - \vec\mu^\top \right) \left(\mathbf \Delta^{-1} \vec{\text {CL}} - \mathbf \Delta^{-1} \vec\mu \right)
\end{aligned}
 \right)}^{- \frac 1 2}
\\ &=
\exp {\left( \begin{aligned}
&\vec p ^\top \mathbf \Sigma^{-1} \vec p -
\vec p ^\top \mathbf \Sigma^{-1} \mathbf X \vec {\text {CL}} -
\vec {\text {CL}} ^ \top \mathbf X ^\top \mathbf \Sigma^{-1}\vec p +
\vec {\text {CL}} ^ \top \mathbf X ^\top \mathbf \Sigma^{-1} \mathbf X \vec{\text {CL}} + \\
& \vec{ \text {CL}} ^\top \mathbf \Delta ^{-1} \vec{\text {CL}} -
 \vec{\text {CL}} ^\top \mathbf \Delta ^{-1} \vec\mu -
\vec \mu ^\top \mathbf \Delta^{-1} \vec{\text {CL}} +
\vec \mu ^\top \mathbf \Delta^{-1} \vec \mu
\end{aligned} \right)}^{- \frac 1 2}
\\ &=
\exp {\left( \begin{aligned}
&\vec{\text {CL}}^\top \left(\mathbf X^\top \mathbf \Sigma^{-1} \mathbf X + \mathbf \Delta^{-1} \right) \vec {\text {CL}} \\
-&\vec{\text {CL}}^\top \left( \mathbf X^\top \mathbf \Sigma^{-1} \vec p + \mathbf \Delta^{-1} \vec \mu \right)
- \left(\vec p ^\top \mathbf \Sigma^{-1} \mathbf X + \vec \mu ^\top\mathbf \Delta^{-1} \right) \vec{\text {CL}} \\
+ & \vec p^\top \mathbf \Sigma^{-1} \vec p + \vec \mu ^\top \mathbf \Delta^{-1} \vec \mu
\end{aligned} \right)}^{- \frac 1 2}
\\ &\propto
\mathcal N \left(
\vec {\text {CL}} \middle | \left(\mathbf X^\top \mathbf \Sigma^{-1} \mathbf X + \mathbf \Delta^{-1} \right) ^{-1} \left( \mathbf X^\top \mathbf \Sigma^{-1} \vec p + \mathbf \Delta^{-1} \vec \mu \right), \left(\mathbf X^\top \mathbf \Sigma^{-1} \mathbf X + \mathbf \Delta^{-1} \right) ^{-1}
\right)
\end{aligned}$$
So, assuming I got my definition of $\mathbf \Delta$ and $\vec \mu$ right, this should be the answer I'm looking for.
