Solution (source code)

= Solution

Condition throughout on the training inputs and $x_*$. Interpret $\theta$ here as the fixed <covariance function> and noise hyperparameters. Assume observation errors $\epsilon_i$ are independent standard <normal random variables>, independent of the <Gaussian process> $f$, and write
$$
K_{ij}=k(x_i,x_j),\quad k_*=(k(x_1,x_*),\ldots,k(x_n,x_*))^T,\quad k_{**}=k(x_*,x_*),\quad C=K+\sigma^2I_n.
$$
The independent-error assumption is the standard <Gaussian process regression> interpretation of the observation model; the single-observation marginal noise law alone would not specify cross-observation noise correlations. We take $\sigma^2>0$ so that $C$ is a <positive-definite matrix>. More generally the inverse formulas below apply whenever the indicated noisy <covariance matrices> are invertible.

The <Gaussian process> prior makes $(f_*,\mathbf f)$ a <Gaussian random vector>. Adding the independent <Gaussian noise> vector to $\mathbf f$ preserves a <multivariate normal distribution>. Its <mean> is zero, its observation <covariance matrix> is $C$, and $\operatorname{Cov}(f_*,y_i)=k(x_*,x_i)$ because the noise is independent of $f_*$. Therefore
$$
\boxed{\begin{pmatrix}f_*\\\mathbf y\end{pmatrix}\biggm|\theta\sim\mathcal N\left(0,\begin{pmatrix}k_{**}&k_*^T\\k_*&C\end{pmatrix}\right)}.
$$
A degenerate latent <variance> is allowed; the displayed <multivariate normal distribution> then denotes a Gaussian law even if it has no full-dimensional <probability density>.