Why does my train data not fall in confidence interval with scikit-learn Gaussian Process? I have a dataset with 1-dimensional input and output and I am applying Gaussian Process regression to it using scikit-learn:
from sklearn import gaussian_process

X, y = get_data()
gp = gaussian_process.GaussianProcess(nugget=0.1)
gp.fit(X, y)
test_X = np.arange(0, 24*60*60, 15*60)
pred_y, pred_mse = gp.predict(test_X, eval_MSE=True)
confidence_interval = np.sqrt(pred_mse) * 2.33

Plotting the result:
import matplotlib.pyplot as plt
plt.scatter(X, y, label='Train data')
plt.fill_between(test_X,
                 pred_y - confidence_interval,
                 pred_y + confidence_interval,
                 alpha=.9, label='98% confidence interval', color='#B22222')
plt.plot(test_X, test_y, label='Prediction')

Shows the following:

This shows that the predicted mean roughly follows the shape of the data, but I would expect that the majority of the blue data points should fall within the red interval, which I believe to show the 98% confidence interval.
Am I wrong in expecting this? Do I need to set different parameters? Increasing the nugget, say to 10, does not make the interval wider, which I also would have expected.
 A: I am not familiar with gaussian_process in scikit-learn so I cannot comment on tool specific aspects. But some general facts about GP might clarify what is going on. 
As in all models, the conclusions (here confidence intervals) only hold if the assumptions (here a GP with a certain kernel and parameters) are met. So maybe your training data is from a process (very) different from the one specified. If this is indeed the case there is no reason for the confidence intervals to hold. You can very this by comparing with training data which is indeed generated according to the assumptions of the GP.
One possible reason why the confidence intervals do not change if you change the nugget size is that the predicition intervals are for the latent function itself, not for its noisy observations. This can be verified by looking into the documentation of gp.predict.
Many tools determine kernel parameters including the nugget by maximising likelihood. So maybe your initial value for the nugget size was just used as a start value and modified in the fitting process. Again you can check this with the documentation or from available diagnostic information. 
A: As far as I can see, the confidence interval that you are calculating is for the mean, not for the actual values, so the points do not need to fall inside it.
A: This will be an answer for scikit-learn 0.18, but I assume something similar was happening in 0.17.
Even though you can get an impression from docs that hyperparameters are optimized automatically, if you keep a default kernel (1.0 * RBF(1.0)), its hyperparameters (=length_scale) are not optimized and neither is noise hyperparameter (=alpha). And since alpha is by default 1e-10, width of the confidence interval will be zero.
To get hyperparameters optimize automatically you have to specify a kernel, e.g. 1.0 * RBF(length_scale=1) + 1.0 * WhiteKernel(). After fitting GP, calling gp.kernel_ returns optimized hyperparameters, e.g. 0.979**2 * RBF(length_scale=1.79) + 0.339**2 * WhiteKernel(noise_level=0.115).
Here's a full example
from sklearn import gaussian_process
from sklearn.gaussian_process.kernels import RBF, WhiteKernel
import matplotlib.pyplot as plt

kernel = 1.0 * RBF(length_scale=1) + 1.0 * WhiteKernel()

X = np.linspace(0, 2 * np.pi, 100)[:, np.newaxis]
y = np.sin(X)[:, 0] + 0.1 * np.random.randn(100)
gp = gaussian_process.GaussianProcessRegressor(kernel=kernel)
gp.fit(X, y)
pred_y, sigma = gp.predict(X, return_std=True)
confidence_interval = sigma * 1.96

plt.scatter(X, y, label='Train data')
plt.fill_between(X[:, 0],
                 pred_y - confidence_interval,
                 pred_y + confidence_interval,
                 alpha=.9, label='98% confidence interval', color='#B22222')
plt.gca().set_title(gp.kernel_)
plt.plot(X, pred_y, label='Prediction')


