value of $Y(\bm{a})$ given the training data. It is a standard result [@ohagan_curve_1978] that $$\begin{aligned}
  \hat{r}(\bm{a})
  &= r(\bm{a}) + \bm{k}^\mathrm{t}(\bm{a})\bm{K}^{- 1} (\bm{y} - \bm{r}), \text{ and}
  \label{eq:gpe_mean} \\
  \hat{\sigma}^2(\bm{a})
  &= k(\bm{a}, \bm{a}) - \bm{k}^\mathrm{t}(\bm{a})\bm{K}^{- 1} \bm{k}(\bm{a}) \label{eq:gpe_var}\end{aligned}$$ where $[\bm{k}(\bm{a})]_i = k(\bm{a}_i, \bm{a})$. These two equations are the principal results of GPE. The estimator, $\hat{r}$, is the sum of two terms. The first is the regression function, and the second a smoothing of the residuals, $(\bm{y} - \bm{r})$. In fact this smoothing is a weighted sum of the residuals, where the weights are the elements of the vector $\bm{k}^\mathrm{t}(\bm{a})\bm{K}^{- 1}$.

Optimizing the emulator {#sec:optimizing_hyperparams}
-----------------------

To evaluate equations \[eq:gpe\_mean\] and \[eq:gpe\_var\], we must assume a mean function, $r$, and a covariance function, $k$. We will choose these from a model, i.e. we will assume that $r \in \{r_{\bm{\mu}} \,|\, \bm{\mu} \in \bm{M}\}$ and $k \in \{k_{\bm{\nu}} \,|\, \bm{\nu} \in \bm{N}\}$, where we call $\bm{M}$ and $\bm{N}$ sets of *hyperparameters*. We optimize our choice of hyperparameters using their maximum likelihood. By equation \[eq:sample\_distrib\] the PDF is $$\begin{aligned}
    \begin{split}
    &f_{\bm{Y}}(\bm{y}; \bm{\theta})\\
    &\quad = \dfrac{1}{\sqrt{(2 \pi)^{N} |\bm{K}|}} \exp \left(-\dfrac{1}{2} (\bm{y} - \bm{r})^{\mathrm{t}} \bm{K}^{-1} (\bm{y} - \bm{r}) \right)
    \end{split}\end{aligned}$$ where the hyperparameter vector $\bm{\theta} := 