# 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.

• @DeltaIV Perhaps I should clarify the first part of the paper's title: is GPR not simply heteroscedastic by adding an input-dependent diagonal to the covariance matrix? Figure 1 appears to make sense; the authors model a given function, $y = (x+w)^2$ by assuming the random variable, $w$, is normally distributed in this instance, then plug that right into the formula for $y$. I suppose I'm just struggling to relate how this is GPR as defined above (i.e. with $y$ following a prescribed mean and variance, it would mean the standard deviation for a predicted $y$ must be symmetric). – Mathews24 Jan 28 at 9:23
• 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 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.

• I believe that's the crux of my question: why does integrating $\ref{a}$ over $p(\theta\vert x,y)$ result in a normal distribution, yet when I introduce a latent variable, $w$, why does integrating over $p(\omega,\theta\vert x,y)$ result in a possibly non-Gaussian distribution? I'm just not quite seeing the difference in the outputs of the two cases mathematically. – Mathews24 Jan 28 at 16:14
• What does it intuitively mean for the data generating process to have additional terms in the argument of the exponential of (11)? In general, the covariance is simply relating the smoothness between points at $x_i$ and $x_j$ based upon a length scale parameter. Can this not handle fitting of arbitrary data? Perhaps to be more precise: what data generating processes could not be modeled with the above assumptions? – Mathews24 Jan 28 at 16:16
• That clears it up tremendously—it's only normally distributed when taking the MLE! But I am wondering then, why is $\theta$ not sufficient for deriving non-Gaussian uncertainty? Why was a latent variable, $w$, needed and why was $\theta$ not enough? – Mathews24 Jan 28 at 16:49