Solution
= Solution
Let $K_\theta$ have entries $k_\theta(t_i,t_j)$ and let
$$
C_\theta=K_\theta+\operatorname{diag}(\sigma_1^2,\ldots,\sigma_N^2).
$$
Independent Gaussian measurement errors give the <multivariate normal density>
$$
\boxed{
p(y\mid t;\mu,\theta)
=(2\pi)^{-N/2}|C_\theta|^{-1/2}
\exp\left[-\frac12(y-\mu\mathbf1)^TC_\theta^{-1}(y-\mu\mathbf1)\right]}.
$$