optimal sequential sampling in gaussian process models Let's say we have a one dimensional dataset of 24 points along with their responses. I am reserving three boundary points for testing (i=1,23,24) and i am fitting a Gaussian process model based on a squared exponential covariance function using GP_fit in R (optimizing the parameters by minimizing the deviance). 
What i am trying to do is to fit the model based on an initial number of training datapoints (let's say n = 5), evaluate the prediction accuracy, and then select the next point based on an optimal sampling criterion which selects the next training point that minimizes the uncertainty in the prediction. 
The problem is that, logically, increasing your training sample size, should lead to better prediction accuracy, but this is not the case here. 
For instance, Figure 1 shows the Gaussian process model based on 11 points .
Figure 2 shows the Gaussian process model based on 17 training points 
The red points in the Figures are the true values of the test points, the blue line is the GP fit, while the red lines are the error bounds.
Note how the prediction accuracy in the GP model of Figure 1 (based on a lesser number of training points) is better. 
Quantitatively, RMSE of GPmodel1 = 0.3405, while RMSE of GPmodel2 = 0.5348.
Is this due to overfitting? because i notice that the GP model fit in both figures goes down significantly after the boundary training points, failing to predict that the original relationship is going up (convex).
If yes/no, what are possible reasons and what are practical ways to overcome this situation? 
Thanks. 
 A: Judging from the pictures your prior mean function seems to be a problem. Far away from data the mean of your posterior/fitted GP will be identical to the prior mean. It looks as if you assumed a prior mean which is a negative constant. This is slightly unusual as the default is often zero-mean. Why did you do that? 
Anyway, as you observe yourself your response seems to be sloping upward which would violate the constant negative mean assumption. Note that your error bounds are only valid, if your model, including the prior mean assumption, is correct.
To remedy the situation you need to choose a different prior mean. You have basically two options: A deterministic function or something stochastic which you include in the calibration.
Choosing a different deterministic mean function is really easy. You just deduct the mean function values from your response and then fit a zero-mean GP to the residuals. 
Having a stochastic model for the mean-function is a possibility if your problem is more complex. You need to find out what your R package supports or code something yourself. The basic idea is described well in Chapter 2.7 of the Rasmussen/Williams book. You parametrise the mean function and assume a prior distribution for those parameters, which can then be treated like any other hyperparameters. For example you can include them in the maximum likelihood optimisation. In equation (2.45) of Rasmussen/Williams you can find the likelihood.
That said, with 17 datapoints you will not have much leeway for additional parameters without overfitting. 
