Variational Inference I have been looking at variational inference, however I came across a problem I have a hard time to understand how to start with it.
Let 
$$\mu_n \sim N(\mu', \sigma'^2) \\ \tau_k \sim N(\tau', \lambda'^2)$$ 
where $\mu', \sigma', \tau', \lambda'$ are given parameters. Now let $$X_n \sim N(\mu_n, \sigma_n^2) \\ Y_k \sim N(\tau_k, \lambda_k^2)$$ where $\lambda_k, \sigma_n$ are given, however the mean $\mu_n$ and $\tau_k$ are stochastic variables (according to the distribution above).
We only observe the sum of the outcomes, i.e. $Z_n^k = Y_k + X_n$, so $Y_k$ and $X_n$ are hidden. Furthermore we observe some data for $n=1:N$, $k=1:K$. Now I want to approximate $p(\mu_{1:N}, \tau_{1:K} | Z_{1:N}^{1:K})$.
Now it want to express the joint probability. However now I get confused which variables should be included. According to the Wikipedia page (section "In practice") they seem to only include the hidden variables of interest and the observed data, i.e. the joint becomes:
$$p(\mu_{1:N}, \tau_{1:K}, Z_{1:N}^{1:K})$$
However, later on (in "A more complex example"), the seem to include all the stochastic variables and the observed data, i.e:
$$p(\mu_{1:N}, \tau_{1:K}, Z_{1:N}^{1:K}, X_{1:N}, Y_{1:K})$$
I was wondering which joint probability is correct to use in this model?
 A: This answer may start a bit too far 'at the beginning', but you can see from which point onwards I start talking about things that are new to you. 
In Bayesian statistics $\mu_n$ and $\tau_k$ are called priors for $X_n$ and $Y_k$ respectively. The priors themselves have priors too, they are called hyper-priors (here $\sigma'$, $\mu'$, $\tau'$ and $\lambda'$). All hyper-priors are constants here. According to the problem statement $\sigma_n$ and $\lambda_k$ may be treated as deterministic as well.
I would say that the joint probability is the probability of all random variables, conditioned on all non-random (hyper)priors. In cases such as these, I like to use my newly learned trick: Graphical models. They allow to visually express the conditional dependencies of probability distributions. Also, the joint probability can be read of relatively easily. The graphical model for your problem looks like this:

The circled entities are random variables, whereas the others are deterministic. With the help of the graphical model, the joint probability is:
\begin{eqnarray}
P(&& X_{1:N}, \mu_{1:N}, Y_{1:K}, \tau_{1:K} | \sigma_{1:N}, \sigma', \mu', \tau', \lambda', \lambda_{1:K}) \\ 
= &&\prod_{n=1}^{N} \underbrace{P(X_n = x_n | \sigma_n, \mu_n)}_{\sim N(\mu_n, \sigma_n)}\underbrace{P(\mu_n | \sigma', \mu')}_{\sim N(\mu', \sigma')} \cdot \prod_{k=1}^{K} \underbrace{P(Y_k = y_k | \tau_k, \lambda_k)}_{\sim N(\tau_k, \lambda_k)}\underbrace{P(\tau_k | \tau', \lambda')}_{\sim N(\tau', \lambda')}\\
= &&P(Z_{n,k} = z_{n,k} | \sigma_{1:N}, \mu', \sigma', \lambda_{1:K}, \tau', \lambda')
\end{eqnarray}
The last equality follows, by the fact that when $X_n$ and $Y_k$ are known, $Z_{n,k}$ is not a random variable anymore. Therefore, we can say that the joint probability here is also the probability distribution of $Z_{n,k}$.
As far as I know (see here for details), the product of Gaussian variables does not result in another Gaussian or something else that's easy to work with (i.e. they are not conjugate priors). Therefore, I think this expression is the most appropriate representation of the joint probability.
