I am trying to use Gaussian Process Regression to model data that has a clear upward rising trend. Here is a basic plot of the data -

enter image description here

The most well-known example of the flexibility of GPs in such problems is shown in the modeling of Mauna Loa CO2 data, examples of which are shown both on the sklearn page and the book "Gaussian processes for machine learning" by Rasmussen. Here is a figure of the fitted model (and predictions) from the sklearn page

sklearn CO2 model (with an upward trend)

However, despite trying various different types of kernels and hyperparameters, I find it hard to model the trend, especially when the GP encounters data from an unseen region.

For example, this code snippet produces the following image -

gp = GaussianProcessRegressor(kernel=kernel,n_restarts_optimizer=100,alpha=1)
testPoints = np.linspace(0,120,120)
yPred,sigma = gp.predict(testPoints.reshape(-1,1),return_std=True)
print("Log-marginal-likelihood: %.3f"% gp.log_marginal_likelihood_value_)
plt.plot(x, y, 'r:')
plt.plot(x, y, 'r.', markersize=10, label=u'Observations')
plt.plot(testPoints, yPred, 'b-', label=u'Prediction')
plt.fill(np.concatenate([testPoints, testPoints[::-1]]),
         np.concatenate([yPred - 1.9600 * sigma,
                        (yPred + 1.9600 * sigma)[::-1]]),
         alpha=.5, fc='b', ec='None', label='95% confidence interval')
plt.legend(loc='upper left')

enter image description here

In order to model the obvious trend outside the 'seen' region in data, one possible solution is illustrated in a Cross Validated post here, which illustrates the benefit of regression-kriging --- essentially learning a polynomial regression model to capture the trend and using a GP model for the residuals. This code snippet for the data shown above yields the following fit (for the regression model) -

from sklearn.preprocessing import PolynomialFeatures
from sklearn import linear_model

poly = PolynomialFeatures(degree=2)
x_ = poly.fit_transform(x)
lm = linear_model.LinearRegression()
yPred = lm.predict(x_)

enter image description here

The following is the result for the combined regression model and a gp fit to the residuals

enter image description here

Here are my questions -

  1. Is it in practice, extremely difficult to model trends with GPs, especially when it entails prediction on a region of unseen data?

  2. In such a case, is regression-kriging the most used method?

  3. How does one calculate the confidence for such a model?


1 Answer 1


I extracted data from the scatter plot and ran it through the "function finder" on my open source Python zunzun.com curve fitting web site. It fit the data to 291 non-linear and 16 linear equations with 3 or less parameters, using a genetic algorithm to estimate initial parameters for the non-linear equations. The idea was to find the simplest smooth curve that gives a reasonably good fit, as that should both interpolate and extrapolate well.

Here is the top result, an inverted offset exponential equation "y = a * exp(b/(x+c))" with parameters a = 1.4271049053537125E+04, b = -1.3705804811393989E+02, and c = 1.2110000000070614E-01 yielding RMSE = 150.2 and R-squared = 0.984. I have extended the plot scaling to the right in order to display potential extrapolation.



Your Answer

By clicking “Post Your Answer”, you agree to our terms of service and acknowledge you have read our privacy policy.

Not the answer you're looking for? Browse other questions tagged or ask your own question.