Fitting a Gaussian Process

Interpolate noisy data with a Gaussian process that returns not just predictions but calibrated uncertainty.

The idea

A Gaussian process (GP) places a distribution over functions. Any finite set of points is jointly Gaussian with covariance given by a kernel. Conditioning on observed data yields a predictive mean and variance at every new input - regression with built-in error bars.

The kernel

The squared-exponential kernel k(x,x') = s^2 exp(-(x-x')^2/(2 l^2)) encodes the belief that nearby inputs have similar outputs. The length scale l sets smoothness; s^2 sets amplitude. Adding observation noise sigma^2 on the diagonal accounts for measurement error.

import numpy as np
def kern(a,b,l=1.0,s=1.0):
    return s*s*np.exp(-(a[:,None]-b[None,:])**2/(2*l*l))
Xtr=np.array([-3,-1,0,2,4.]); ytr=np.sin(Xtr)
Xte=np.linspace(-5,6,100); noise=1e-4
K=kern(Xtr,Xtr)+noise*np.eye(len(Xtr))
Ks=kern(Xtr,Xte); Kss=kern(Xte,Xte)
Kinv=np.linalg.inv(K)
mu=Ks.T@Kinv@ytr
var=np.diag(Kss-Ks.T@Kinv@Ks)
print('pred at 1.0:',round(float(mu[np.argmin(abs(Xte-1))]),3))
print('max std:',round(float(np.sqrt(var.max())),3))

The uncertainty

The predictive variance is small near observed points and grows in the gaps and beyond the data - exactly the honest behaviour you want. This makes GPs ideal for Bayesian optimization, where you decide where to sample next by trading off predicted value against uncertainty (exploration versus exploitation).

Costs and choices

Training requires inverting an n-by-n kernel matrix, order n^3, so plain GPs suit small-to-medium datasets; sparse and inducing-point approximations extend them further. The hyperparameters l, s, and noise are usually fit by maximizing the marginal likelihood, which automatically balances data fit against smoothness.