its task. Here we restrict ourselves to using the squared-exponential covariance function (equation \[eq:cov\_SE\]), and hence consider the second approach.

We have assumed that our data is drawn from a Gaussian random process with covariance $k_{\mathrm{SE}}(\bm{a}, \bm{a}') + \sigma_\epsilon^{2}(\bm{a})$ (equation \[eq:gpe\_assumption\]). Crucially, the squared exponential covariance function is a function of the difference of its arguments only and the properties of the random process are therefore independent of the absolute values its arguments. In particular, the variance and autocovariance of the random process are constant. The fact that the variance is constant means that noise-free data are identically normal. If we do not observe this distribution in our training data we may transform it to ensure that we do. Such a transformation is said to be *variance stabilizing*.

One such variance-stabilizing transformation is the *Box-Cox transformation*. Let $\{y(\bm{a}_i)\}_{i = 1}^N$ be a set of observations. Then the Box-Cox transformation of the observations [@box_analysis_1964] is the function $g$ such that $$\begin{aligned}
  g(y(\bm{a}_i); \lambda_1, \lambda_2) = \dfrac{(y(\bm{a}_i) + \lambda_2)^{\lambda_1} -
    1}{\lambda_1}\end{aligned}$$ for some real $\lambda_1$, $\lambda_2$ such that $\lambda_{2} > - y(\bm{a}_{i})$ for all $i$. Note that this expression is just a scaled power law, with the scaling chosen such that $\lim_{\lambda_1 