# Using Gaussian process regression with non Gaussian data

I have a question about practical implementation and interpretation of the Gaussian process regression model given by Rasmussen & Williams.

The regression problem is defined as follows:

Let $$x_i\in \mathbb R^{5}$$ be an input vector and $$y_i\in \mathbb R$$ be its corresponding target. The set of $$M$$ inputs are arranged into a matrix $$X=[x_1,\ldots,x_M]^⊤$$ and their corresponding targets are stored in a matrix $$Y=[y_1,\ldots, y_M]^\top.$$

I wish to train a GPR model $$G={X,Y,\theta}$$ using the squared exponential function:

$$k(x_i,x_j)=\alpha^2\exp\left(\frac{−1}{2\beta^2(x_i−x_j)^2}\right)+\gamma^2\delta_{ij},$$

where $$\delta_{ij}$$ equals 1 if $$i=j$$ and 0 otherwise. The hyperparameters are $$\theta=(\alpha,\beta,\gamma)$$ with γγ being the assumed noise level in the training data and $$\beta$$ is the length-scale.

To train the model, I need to minimise the negative log marginal likelihood with respect to the hyperparameters:

$$−\log p(Y∣X,\theta)=\frac{1}{2}\operatorname{tr}\left(Y^\top K^{−1}Y\right)+\frac{1}{2}\log∣K∣+c,$$ where $$c$$ is a constant and the matrix $$K$$ is a function of the hyperparameters

covfunc = @covSEiso;
likfunc = @likGauss;
sn = 0.1;
hyp.lik = log(sn);
hyp2.cov = [0;0];
hyp2.lik = log(0.1);
hyp2 = minimize(hyp2, @gp, -100, @infExact, [], covfunc, likfunc, X1, Y1(:,n));
exp(hyp2.lik)
nlml2 = gp(hyp2, @infExact, [], covfunc, likfunc, X1, Y1(:, n));
[m s2] = gp(hyp2, @infExact, [], covfunc, likfunc, X1, Y1(:, n), XT);
YT(:, n) = m;


X1,Y1 are the training inputs/targets

XT are the test inputs, YT are the target predictions.

My question is about the use of [m,s2] as predictions for Y (the second to last line of code). If the original sample Y are not normally distributed is the predictive distribution given by the values in [m,s2] stil valid? For example is it valid to say that for a test value $$x_{*i}$$ the distribution of the prediction for the target $$y_{*i}$$ is $$N(m,s2)$$, or am I misinterpreting the implementation of this model?

• In Gaussian process regression (GPR) we assume that our function is a realization of Gaussian process. Almost any function is close to some Gaussian process realization. In particular, GPR can fit the training data perfectly i.e. pass through the all sample points. I have encountered no discussion regarding testing of Gaussian nature of data while using GPR. Sep 29, 2016 at 11:43
• The assumption in GPR is that the noise is Gaussian. The sample Y will not typically be, but the residuals after you've applied your model should be if the model is well-specified. This is illustrated in the figure below, somewhat. There are cases for which the noise is truly non-gaussian so that even with a well-specified model the residuals will be non gaussian or heteroskedastic. In these cases, something like Warped Gaussian Processes by Snelson/Ghahramani can be used. There is much literature on the topic, they provide a few good references towards the end of the paper. Feb 4, 2019 at 18:17
• @DavidKozak are there python or R packages that can be used to implement warped gaussian processes? Aug 31, 2021 at 18:44
• @AngusCampbell I believe that GPy in Python does. Sep 1, 2021 at 20:20

For a gaussian process, the evaluation of the gaussian processes of an input value $(x+1)$ is the output $y(x+1)$ obtained by
$$\mathord{\buildrel{\lower3pt\hbox{\scriptscriptstyle\frown}} \over y} = {K^y}(\mathord{\buildrel{\lower3pt\hbox{\scriptscriptstyle\frown}} \over x} ,x)*{\left( {{K^Y}} \right)^{ - 1}}*y$$
in your code, this is expressed by $m$ ,the mean value is the value of higher probability. For each value that you want to evaluate, you need to calculate $${K^y}(\mathord{\buildrel{\lower3pt\hbox{\scriptscriptstyle\frown}} \over x} ,x)$$ and then you apply the equation for the mean value to obtain your output value.