# How to prevent heteroskedastic models from overfitting?

I'm trying to fit neuroscience data using a Gaussian Process, but noticed that it behaves poisson-like (var = mean).

Since classic GP models assume iid noise, I figured I could get a better fit by using a heteroskedastic model. For now, I am still assuming a gaussian likelihood, but not sure whether this is valid (what are the customary or usual tests that people use to check whether the data is gaussian-like?)

But if we go along this assumption, I made some simple tests to see how heteroskedastic GP models would perform, and here are the results on a test function (with homoskedastic normal errors):

The heteroskedastic model has the best negative log marginal likelihood (nlml) but is clearly not fitting the data properly. My questions are then the following:

1. Aren't bayesian models self-regularizing and hence don't overfit? Mackay explained it with the following picture, which I got from p. 110 of the Gaussian Process for Machine Learning textbook.

Perhaps this is only true if we do the full bayesian model, whereas here we are doing MAP inference on the hyperparameters?

1. How should I prevent this overfitting?

There are many methods possible. One is to change the likelihood. Another is different models for predicting the mean of the variance term: from simple linear model, to ARCH models, to full-out GP model (which I believe is what the GPy code I am using is doing).

Here is the code if anyone wants to play with it:

import numpy as np
import pylab as pb
import GPy
import matplotlib.pyplot as plt

np.random.seed(0)

def f(X):
return  10. + .1*X + 2*np.sin(X)/X

X = np.random.uniform(-10,20, 50)
X = X[~np.logical_and(X>-2,X<3)] #Remove points between -2 and 3 (just for illustration)
X = np.hstack([np.random.uniform(-1,1,1),X]) #Prepend a point between
#-1 and 1  (just for illustration)
#X = np.linspace(-10,20, 50)
error = np.random.normal(0,.2,X.size)
# error = np.zeros(X.size)
# for i in range(X.size):
#     error[i] = np.random.normal(0,abs(f(X[i])))
Y = f(X) + error

#Heteroscedastic model fixed variance (true value)
kern = GPy.kern.MLP(1) + GPy.kern.Bias(1)
m = GPy.models.GPHeteroscedasticRegression(X[:,None],Y[:,None],kern)
m['.*het_Gauss.variance'] = abs(error)[:,None] #Set the noise parameters to the error in Y
m.het_Gauss.variance.fix() #We can fix the noise term, since we already know it
m.optimize()

#Heteroscedastic model fixed variance (0.05)
kern = GPy.kern.MLP(1) + GPy.kern.Bias(1)
m1 = GPy.models.GPHeteroscedasticRegression(X[:,None],Y[:,None],kern)
m1.het_Gauss.variance = .05
m1.het_Gauss.variance.fix()
m1.optimize()

#Heteroscedastic model (learn variance)
kern = GPy.kern.MLP(1) + GPy.kern.Bias(1)
m2 = GPy.models.GPHeteroscedasticRegression(X[:,None],Y[:,None],kern)
m2.optimize()

# Homoscedastic model (learn variance)
kern = GPy.kern.MLP(1) + GPy.kern.Bias(1)
m3 = GPy.models.GPRegression(X[:,None],Y[:,None],kern)
m3.optimize()

print(m)
print(m1)
print(m2)
print(m3)

fig,axs = plt.subplots(4,1, sharex=True, sharey=True)
m.plot(ax=axs[1])
axs[1].set_title("Heteroscedastic fixed variance (= sampled error): nlml = {}".format(m.objective_function()))

m2.plot(ax=axs[2])
axs[2].set_title('Heteroscedastic model (learn variance): nlml = {}'.format(m2.objective_function()))

m3.plot(ax=axs[0])
axs[0].set_title('homoscedastic model (learn variance): nlml = {}'.format(m3.objective_function()))

axs[3].plot(np.linspace(-10,20),f(np.linspace(-10,20)),'r-')
axs[3].plot(X,Y,'kx',mew=1.5)
axs[3].set_title('real function + observations')
plt.show()