GPFlow: how to account for uncertainties from mean model

Question

In GPFlow one can add a fitted mean function to the GP regression. When doing this as in the basic example, the result is, that there will be no uncertainties due to the uncertainty in the fit of the mean. E.g. in the example below the error bars don't grow outside the range of available data, as the slope of the linear mean remains fixed at its optimized value. Is there a way to account for these uncertainties, such that the error bands grow when extrapolating?

(The question was originally stated in an issue report but moved here to be more accessible)

import numpy as np
import matplotlib.pyplot as plt
import gpflow
from gpflow.utilities import print_summary

def f(x):
    return np.sin(3*x) + x

xtrain = np.linspace(0, 3, 50).reshape([-1, 1])
ytrain = f(xtrain) + 0.5*(np.random.randn(len(xtrain)).reshape([-1, 1]) - 0.5)

k = gpflow.kernels.SquaredExponential()
meanf = gpflow.mean_functions.Linear()
m = gpflow.models.GPR(data=(xtrain, ytrain), kernel=k, mean_function=meanf)
opt = gpflow.optimizers.Scipy()

def objective_closure():
    return - m.log_marginal_likelihood()

opt_logs = opt.minimize(objective_closure,
                        m.trainable_variables,
                        options=dict(maxiter=100))

print_summary(m)


xpl = np.linspace(-5, 10, 100).reshape(100, 1)
mean, var = m.predict_f(xpl)

plt.figure(figsize=(12, 6))
plt.plot(xtrain, ytrain, 'x')
plt.plot(xpl, mean, 'C0', lw=2)
plt.fill_between(xpl[:, 0],
                 mean[:, 0] - 1.96 * np.sqrt(var[:,0]),
                 mean[:, 0] + 1.96 * np.sqrt(var[:,0]),
                 color='C0', alpha=0.2)

Answer 1

Most of GPflow's models only optimise for the MAP estimate of the hyperparameters of the kernel, mean function and likelihood. The models do not account for uncertainty on these hyperparameters during training or prediction. While this could be limiting for certain problems, we often find that this is a sensible compromise between computational complexity and uncertainty quantification.

That being said, in your specific case (i.e. a linear mean function) we can account for uncertainty in the linear trend of the data by specifying a linear kernel function, rather than a linear mean function.

Using your snippet with this model specification:

k = gpflow.kernels.SquaredExponential() + gpflow.kernels.Linear()
meanf = gpflow.mean_functions.Zero()  
m = gpflow.models.GPR(data=(xtrain, ytrain), kernel=k, mean_function=meanf)

Gives the following fit, with error bars that grow outside the data range: