Gaussian Process: why is my data only explained by noise? I'm trying to explain a very simple 10-pt dataset by a Gaussian Process (part of a larger bayesian optimization framework) and I don't understand why it is only being explained by noise.
Here is the data, on a 2x5 grid.

I am able to perfectly fit the data if only using an RBF kernel, which gets optimized to 

6.54**2 * RBF(length_scale=[0.01, 0.01])

But if I add a noise WhiteKernel, then the data gets explained by the noise

4.85**2 * RBF(length_scale=[146, 6.6]) + WhiteKernel(noise_level=17.2)

and the mean prediction looks like:

I need to add noise because the real data I am working with is noisy (EMG responses to brain stimulation are not deterministic). However I noticed it was mostly explained by noise, so tried to see what the GP would learn with the mean (deterministic) function of the data, and it is again only explained by noise.
I do know about the 5th chapter of Rasmussen and Williams' Gaussian Process Textbook, and the fact that the marginal likelihood can be multimodel, since some functions can equally well be explained by noise or quickly-varying functions, but in this case only the noise version seems to happen. Here is the code for reference.
import matplotlib.pyplot as plt
from mpl_toolkits.mplot3d import Axes3D

x = np.arange(0,2)
y = np.arange(0,5)
xx,yy = np.meshgrid(x,y)
xr,yr = xx.ravel(), yy.ravel()
xtrain = np.c_[xr,yr]

ytrain = np.array([[ 0.87611],
                   [ 3.58377],
                   [ 1.32999],
                   [ 7.63373],
                   [ 1.09371],
                   [ 9.18058],
                   [13.91528],
                   [ 3.21335],
                   [ 3.26006],
                   [ 7.30752]])

fig = plt.figure()
ax = fig.gca(projection='3d')
surf = ax.scatter(xx, yy, ytrain.reshape((5,2)), antialiased=False)

from sklearn.gaussian_process import GaussianProcessRegressor
from sklearn.gaussian_process.kernels import RBF, WhiteKernel

kernel = 1.0 * RBF(length_scale=[1.,1.],
                   length_scale_bounds=[(1e-2, 1e3), (1e-2, 1e3)]) \
    + WhiteKernel(noise_level=1, noise_level_bounds=(1e-10, 1e+2))
gp = GaussianProcessRegressor(kernel=kernel,
                              alpha=1e-10,
                              n_restarts_optimizer=5).fit(xtrain, ytrain)
print(gp.kernel_)

fig = plt.figure()
ax = fig.gca(projection='3d')
surf = ax.scatter(xx, yy, gp.predict(xtrain).reshape((5,2)), antialiased=False)
plt.show()

 A: I've played around with this a little bit, and it looks like you have two local minima: one explains almost everything with noise, and the other explains almost everything by high-frequency oscillations. This is easier to see if you do a couple of things: 


*

*Normalize your training data to have mean 0 and variance 1. Because your GP has a zero-mean prior, shifting your data to 0 mean helps. Normalizing the variance to 1 makes the noise variance more interpretable. You can do this by setting


gp = GaussianProcessRegressor(kernel=kernel,
                              alpha=1e-10,
                              n_restarts_optimizer=5,
                              normalize_y=True).fit(xtrain, ytrain)



*Change the bounds on your hyperparameters to reflect this. You'll need to set the lower bounds on the lengthscale much lower (I used 1e-10), and I'd upper bound the noise at 1 (it doesn't make sense for the noise variance to be greater than the data variance). I'd also do 100 restarts instead of 5. 


Even after doing this, the hyperparameters are pretty unstable: I get different results running it multiple times, with a cluster where the noise is 0.3 - 0.6 and the lengthscales are 1e-10, and another where the length scales are more moderate but the noise is 1e-6 - 1e-9. If you want to a more complete picture of what's going on, you probably need to set proper priors over the hyperparameters and use Hamiltonian Monte Carlo to estimate the posterior. This is a pretty good guide. 
