Including error dependent on output in Gaussian Process Regression

Question

I have a set of experimental data that I am trying to fit using Gaussian process regression (GPR) using Python's sklearn package. The only problem is that my data has an experimental standard deviation equal to 20% of its output. Can this dependency be handled in GPR?

In the code, the essential aspects are:

total_noise = 0.2*((numpy.abs(y))) #standard deviation

# Instantiate a Gaussian Process model
gp = GaussianProcessRegressor(kernel=kernel, alpha = total_noise**2.,
                                  n_restarts_optimizer=100)

If my kernel is simply a radial basis function, then this results in:

$$k(x_i, x_j) = RBF(x_i, x_j) + 0.2 *y_i * \mathcal I (x_i * x_j)$$

But the problem appears to be that alpha here is dependent on y (my experimental output). Correct me if I am mistaken, but I believe this is okay when training. But I ultimately want both the computed mean and standard deviation from GPR, and I do not believe sklearn's GPR is accounting for this dependency on the output in alpha since I would expect the predicted standard deviation, $\sigma^*$, to be smaller where the predicted output, $y^*$, is smaller. But this is not the case. My covariance itself appears to depend on $y$ because, in the model, $y$ exists on both sides of the sampling statement:

$$y \sim \mathcal N(0, \Sigma(y))$$

and I do not believe this is handled consistently. Is there a way to account for this dependency (i.e. my original data has a standard deviation of 20% in the y-axis) using sklearn's GPR?

Answer 1

If the error bars represent observation uncertainty, then essentially, you're trying to model a latent function $\bar Y$ which is a function of $X$, which is corrupted by noise:

$$Y = \bar Y + \bar Yk\epsilon$$

... where $k\epsilon \sim \mathcal N(0, k^2 \sigma^2)$.

Notice that this is quite a complicated model, i.e. if one models $\bar Y$ using a GP, the distribution of $\bar Y k \epsilon$ is no longer multivariate normal, hence Y isn't a GP. You could easily model this using a probabilistic programming tool such as Stan or PyMC, but given that $Y$ isn't a GP, I do not think that it's possible to fit the model above in scikit-learn.

I'm not entirely sure if adding an alpha term in scikit-learn is a good way to approximate the solution (I think that it is, as the variance corresponding to observed points is higher, but perhaps this is something others can weigh on).

Note that the predicted error bars of the GP don't necessarily need to look like the error bars you drew from your data - you've got lots of data, so even though the variance of the initial data points is high, there are many of them, which gives the GP enough information as to where the function will lie, roughly. The model as it's written above does not correspond well to the data, as even though the error bars are smaller along the domain, notice that they're all over the place. Really, the variance of the data hasn't really changed all that much over the domain, so there's no reason to expect the GP to have smaller error bars further along the line.

(Look at the point near (20, 60) - I do not see any reasonable line passing through all the points while exactly maintaining the error bar 'profile' as in that graph. There's more to the errors than the model above, so a GP is a very reasonable approximation to the mean function of that graph.)