Predictive posterior distribution with multivariate normal distribution Suppose I have a multivariate normal ${\bf{Y}}|{\bf{\theta}} \sim {\bf{MVN}}(X {\bf{\beta}}, \sigma^{2}H(\phi))$ where ${\bf{Y}}$ is a set of observations ${\bf{Y}} = \{y({\bf{s}}_{1}),y({\bf{s}}_{2}),... ,y({\bf{s}}_{n})\}$ and $H(\phi)$ is a covariance matrix $H(\phi) = \rho({\bf{s}}_{i} - {\bf{s}}_{j};\phi)$. Now I want to predict a value $y_{0}$, so as it is usual, the predictive posterior distribution is
$$p(y_{0}|{\bf{y}}, X, x_{0}) = \int p(y_{0}|{\bf{y}},{\bf{\theta}}, x_0)p({\bf{\theta}}|{\bf{y}}, X)d\theta$$
Usually, this is described as a product of $\text{Likelihood} \times \text{Posterior}$. However, in this case, the likelihood was initially expressed as a multivariate normal but I just want to predict the value $y({\bf{s}})$ at some point ${\bf{s}}_{0}$, so obviously I can't express a single point as a multivariate normal. What should I do? 
As an alternative, I should be able to use a conditional distribution $y_{0} |  \bf{Y}$ from a joint distribution $\left( y_{0}, {\bf{Y}} \right)$ in the following way:
$$\left( \begin{array}{ccc}
Y_{1}  \\
Y_{2}
 \end{array} \right) \sim N\left(\left( \begin{array}{ccc}
\mu_{1}  \\
\mu_{2}
 \end{array} \right), \left( \begin{array}{ccc}
\Omega_{11} & \Omega_{12}  \\
\Omega_{21} & \Omega_{22}
 \end{array} \right)\right)$$
The problem is that this in some way breaks with the idea of using $\text{Likelihood} \times \text{Posterior}$ for the predictive posterior distribution.
Thanks!
 A: It sounds like you are trying to do inference with a Gaussian process.  The correct approach is to use the conditional distribution, as in the second part of your question.  The reason that you cannot use the Likelihood $\times$ Posterior approach is because $\theta$ is an infinite-dimensional object, representing the value of the process at all locations ${\bf s}$.
A: I think I got a very simple answer.
I overlooked the fact that the predictive posterior distribution is derived with a "conditional likelihood":
\begin{array} 
pp(y_{0}|{\bf{y}}) &= \int p(y_{0}, \theta| {\bf{y}})d\theta \\ 
& = \int p(y_{0}|{\theta, \bf{y}})p(\theta| {\bf{y}})d\theta
\end{array}
where I omitted ${\bf{X}}$ and $x_{0}$ for simplicity. The full predictive posterior distribution is what I initially wrote in my question:
$$p(y_{0}|{\bf{y}}, X, x_{0}) = \int p(y_{0}|{\bf{y}},{\bf{\theta}}, x_0)p({\bf{\theta}}|{\bf{y}}, X)d\theta$$
Only when $y_{0}$ is independent of the previous observations ${\bf{y}}$, it is possible to set $p(y_{0}|{\bf{y}},{\bf{\theta}}, x_0) = p(y_{0}|{\bf{\theta}}, x_0)$ to obtain the regular form $$\text{Predictive} = \text{Likelihood} \times \text{Posterior}$$
$$p(y_{0}|{\bf{y}}, X, x_{0}) = \int p(y_{0}|{\bf{\theta}}, x_0)p({\bf{\theta}}|{\bf{y}}, X)d\theta$$
