Using Gaussian process regression with non Gaussian data

Question

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?

Answer 1

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.

You can see a gaussian process as a sum of gaussian distributions over a specific path.