# Deriving non-Gaussian uncertainties using Gaussian process regression

I came across this work on "Gaussian Process Regression with Heteroscedastic or Non-Gaussian Residuals" which is intriguing, but I believe I am failing to properly interpret its results. Namely, it appears non-Gaussian uncertainties can be derived using Gaussian process regression (GPR), but I am not recognizing how that is possible.

For example, in Figure 2, a plot is produced whereby the residuals are clearly asymmetric and non-Gaussian after running multiple samples ($$n = 19$$ to be specific). But how can this be the case in general (i.e. $$n \rightarrow \infty$$)? Intuitively, how does a latent or unobserved variable permit non-Gaussian uncertainties in GPR? The conditional mean, $$\mu_{y^* \lvert x,y,x^*}$$, and variance, $$\Sigma_{y^* \lvert x,y,x^*}$$, of the prediction is given by (13) and (14), respectively. Is not the fundamental assumption of GPR that:

$$y^* \lvert x,y,x^* \sim N(\mu_{y^* \lvert x,y,x^*},\Sigma_{y^* \lvert x,y,x^*})?$$

Therefore, in Figure 2, how is it possible that the predicted output, $$y^* \lvert x,y,x^*$$, can have residuals on average that are non-Gaussian? In general, can any data set with non-Gaussian errors be modeled using this approach or must the hidden variable, $$w$$, be known?

Any references to GPR codes/numerical packages incorporating non-Gaussian residuals and further literature on the topic would be appreciated as I cannot find anything explicitly beyond this author's own work.

• In the cases in the paper, the expected function for $y$ appears to always be known. But given observed data from an unknown sample with possibly non-Gaussian errors on the observed data, could this paper's method be applied as GPR normally can if the observed errors were Gaussian? – Mathews24 Jan 28 '19 at 9:26

The conditional mean, $$\mu_{y^* \lvert x,y,x^*}$$, and variance, $$\Sigma_{y^* \lvert x,y,x^*}$$, of the prediction is given by (13) and (14), respectively. Is not the fundamental assumption of GPR that:

$$y^* \lvert x,y,x^* \sim N(\mu_{y^* \lvert x,y,x^*},\Sigma_{y^* \lvert x,y,x^*})?$$

The model presented in Wang & Neal is not GPR. It's one of the hundreds of models in use for GPR with heteroskedastic residuals, and in particular it's the GPLV model. The fundamental assumption of GPLV is that

$$y^* \lvert \omega, \theta, x,y,x^* \sim N(\mu_{y^*,\omega, \theta \lvert x,y,x^*},\Sigma_{y^*,\omega, \theta \lvert x,y,x^*}) \tag{a}\label{a}$$

to put it in your terms (to be more precise, the fundamental assumptions are actually eq. (11) of Wang & Neal, and an uncorrelated multivariate standard normal prior on $$\omega_1,\dots,\omega_n$$). When you integrate $$\ref{a}$$ with respect to $$p(\omega,\theta\vert x,y)$$, there's no reason why $$y^*$$ should still have a normal distribution. You can indeed compute mean and variance of $$y^*$$ with the expressions (13) and (14) of Wang & Neal, but there's no reason why the higher order moments should be 0. See the example in Fig.1 of the paper.

You furthermore ask in the comments whether this GPLV model may or may not be applicable to cases where the "expected function" for $$f$$ is unknown. The answer is simple: as long as your data generating process can be modeled in the following hierarchical way

\begin{align} y&\sim f(x) + \epsilon \\\\ f(x)&\sim\mathcal{GP}(\mu(x),k(x,x'; \theta,\omega)),\quad \epsilon \sim\mathcal{N}(0,\sigma^2) \\\\ \theta &\sim ??, \quad \omega \sim \mathcal{N}(0,1) \end{align}

Note that I haven't specified the prior for $$\theta$$, because it's the usual one for a SE kernel. Also, for simplicity I wrote the kernel as $$k(x,x';\theta,\omega)$$, but the GPLV model actually uses a well-defined, specific kernel.

If you know how data are generated for your problem, you can decide if this model is realistic or not. If you don't, then you have a bigger problem than not simply ignoring whether GPLV is the right tool for you or not. You could still try to use the GPLV model, and then apply the vast array of posterior predictive checks developed for Bayesian models in the last 20 years, to verify if the underlying assumptions are reasonable.