the regression function. In Gaussian-process regression we assume that $$\begin{aligned}
    Y \sim {\operatorname{GP}}(r, \sigma^2)
    \label{eq:gpe_assumption}\end{aligned}$$ whereupon $r = {\operatorname{E}}(Y)$, as required. Then $E(\bm{a}) = Y(\bm{a}) -
r(\bm{a})$, i.e. $E$ is a Gaussian process with zero mean. We will assume that its covariance is the sum of two terms, i.e. that $\sigma^{2}(\bm{a},
\bm{a}') = k(\bm{a}, \bm{a}') + \sigma^{2}_{\epsilon}(\bm{a})$, the first term representing the signal and the second term representing the noise (i.e. the error in the dependent variable). Note that the noise is not in general constant, but is a function of $\bm{a}$. In statistical parlance, we would say that the errors are *heteroscedastic* rather than *homoscedastic*.

Suppose we wish to predict the dependent variable $Y(\bm{a})$ for some value $\bm{a} \in \bm{A}$. By hypothesis, the distribution of this random variable is univariate normal, i.e.$$\begin{aligned}
    Y(\bm{a}) \sim N(r(\bm{a}), \sigma^{2}(\bm{a}, \bm{a})),\end{aligned}$$ and the joint distribution of the sample is multivariate normal, i.e. $$\begin{aligned}
    \bm{Y} \sim N(\bm{r}, \bm{K}),
    \label{eq:sample_distrib}\end{aligned}$$ where $\bm{r} := (r(\bm{a}_{1}), r(\bm{a}_{2}), \ldots, r(\bm{a}_{N}))$, and $[\bm{K}]_{ij} = k(\bm{a}_{i}, \bm{a}_{j}) + \sigma_{\epsilon}^{2}(\bm{a}_{i})\delta_{ij}$. The estimator for our regression function is the mean of the conditional random variable $$\begin{aligned}
    Y(\bm{a}) \,|\, (\bm{Y} = \bm{y}) \sim N(\hat{r}, \hat{\sigma}^2)
    \label{eq:gpe_conditional}\end{aligned}$$ where $[\bm{y}]_{i} = y(\bm{a}_{i})$, i.e. the expected 