# Using Gaussian process regression with non Gaussian data

I have a question about practiacal implementation and interpretation of the Gaussian process regression model given by Rasmussen & Williams (http://www.gaussianprocess.org/gpml).

The regression problem is defined as follows:

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

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

$k(x_i,x_j)=α^2exp(\frac{−1}{2β^2(xi−xj)^2})+γ^2δ_{ij}$,

where $δ_{ij}$ equals 1 if $i=j$ and 0 otherwise. The hyperparameters are $θ=(α,β,γ)$ with γγ being the assumed noise level in the training data and $β$ is the length-scale.

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

$−logp(Y∣X,θ)=\frac{1}{2}tr(Y^⊤K^{−1}Y)+\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. – Alexey Zaytsev Sep 29 '16 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. – David Kozak Feb 4 '19 at 18:17

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.