Gaussian Process Regression

  • Author: Johannes Maucher
  • Last Update: 4th December 2013
In [1]:
import numpy as np
from scipy import r_
from matplotlib import pyplot as plt
np.set_printoptions(precision=3,suppress=True)

Definition of training data

In [2]:
xB=np.array([1,2,3,4])
yB=np.array([0.25,0.95,2.3,3.9])

Define the positions at wich the function values shall be predicted

In [3]:
xPred=np.arange(0,7,0.2)

Define the hyperparameters for the mean- and covariance function

In [4]:
c2=0.25 # constant coefficient of quadratic term in prior mean function
ell=1 # horicontal length scale parameter in the squared exponential function
sigmaF2=2 #sigmaF2 is the variance of the multivariate gaussian distribution
sigmaN2=0.005 #sigmaN2 is the variance of the regression noise-term

Definition of mean- and covariance function. Here, the mean function is a quadratic polynomial and the covariance function is the squared exponential.

In [5]:
def priormean(xin):
    return c2*xin**2

def corrFunc(xa,xb):
    return sigmaF2*np.exp(-((xa-xb)**2)/(2.0*ell**2))

Calculate the values mx of the mean function in the range from 0 to 7. These values are just used for plotting. The values mxB are the mean-function values at the training data x-values, i.e. the mean-vector. These values are applied for calculating the prediction.

In [6]:
x=np.arange(0,7,0.1)
mx=priormean(x)
mxB=priormean(xB)

Calculate the covariance matrix by evaluating the covariance function at the training data x-values.

In [7]:
KB=np.zeros((len(xB),len(xB)))
for i in range(len(xB)):
    for j in range(i,len(xB)):
        noise=(sigmaN2 if i==j else 0)
        k=corrFunc(xB[i],xB[j])+noise
        KB[i][j]=k
        KB[j][i]=k        
print '-'*10+' Matrix KB '+'-'*10
print KB.round(decimals=3)
---------- Matrix KB ----------
[[ 2.005  1.213  0.271  0.022]
 [ 1.213  2.005  1.213  0.271]
 [ 0.271  1.213  2.005  1.213]
 [ 0.022  0.271  1.213  2.005]]

Calculate the inverse of the covariance matrix

In [8]:
KBInv=np.linalg.inv(KB)
print '-'*10+' Inverse of Matrix KB '+'-'*10
print KBInv.round(decimals=3)
---------- Inverse of Matrix KB ----------
[[ 0.953 -0.862  0.519 -0.208]
 [-0.862  1.688 -1.219  0.519]
 [ 0.519 -1.219  1.688 -0.862]
 [-0.208  0.519 -0.862  0.953]]

Calculate the covariance matrix $K_*$ between training x-values and prediction x-values

In [9]:
Ks=np.zeros((len(xPred),len(xB)))
for i in range(len(xPred)):
    for j in range(len(xB)):
        k=corrFunc(xPred[i],xB[j])
        Ks[i][j]=k
print '-'*10+' Matrix Ks '+'-'*10
print Ks.round(decimals=3)
---------- Matrix Ks ----------
[[ 1.213  0.271  0.022  0.001]
 [ 1.452  0.396  0.04   0.001]
 [ 1.671  0.556  0.068  0.003]
 [ 1.846  0.751  0.112  0.006]
 [ 1.96   0.974  0.178  0.012]
 [ 2.     1.213  0.271  0.022]
 [ 1.96   1.452  0.396  0.04 ]
 [ 1.846  1.671  0.556  0.068]
 [ 1.671  1.846  0.751  0.112]
 [ 1.452  1.96   0.974  0.178]
 [ 1.213  2.     1.213  0.271]
 [ 0.974  1.96   1.452  0.396]
 [ 0.751  1.846  1.671  0.556]
 [ 0.556  1.671  1.846  0.751]
 [ 0.396  1.452  1.96   0.974]
 [ 0.271  1.213  2.     1.213]
 [ 0.178  0.974  1.96   1.452]
 [ 0.112  0.751  1.846  1.671]
 [ 0.068  0.556  1.671  1.846]
 [ 0.04   0.396  1.452  1.96 ]
 [ 0.022  0.271  1.213  2.   ]
 [ 0.012  0.178  0.974  1.96 ]
 [ 0.006  0.112  0.751  1.846]
 [ 0.003  0.068  0.556  1.671]
 [ 0.001  0.04   0.396  1.452]
 [ 0.001  0.022  0.271  1.213]
 [ 0.     0.012  0.178  0.974]
 [ 0.     0.006  0.112  0.751]
 [ 0.     0.003  0.068  0.556]
 [ 0.     0.001  0.04   0.396]
 [ 0.     0.001  0.022  0.271]
 [ 0.     0.     0.012  0.178]
 [ 0.     0.     0.006  0.112]
 [ 0.     0.     0.003  0.068]
 [ 0.     0.     0.001  0.04 ]]

Calculate the covariance matrix $K_{**}$ between prediction x-values

In [10]:
Kss=np.zeros((len(xPred),len(xPred)))
for i in range(len(xPred)):
    for j in range(i,len(xPred)):
        noise=(sigmaN2 if i==j else 0)
        k=corrFunc(xPred[i],xPred[j])+noise
        Kss[i][j]=k
        Kss[j][i]=k
print '-'*10+' Matrix Kss '+'-'*10
print Kss.round(decimals=3)
---------- Matrix Kss ----------
[[ 2.005  1.96   1.846 ...,  0.     0.     0.   ]
 [ 1.96   2.005  1.96  ...,  0.     0.     0.   ]
 [ 1.846  1.96   2.005 ...,  0.     0.     0.   ]
 ..., 
 [ 0.     0.     0.    ...,  2.005  1.96   1.846]
 [ 0.     0.     0.    ...,  1.96   2.005  1.96 ]
 [ 0.     0.     0.    ...,  1.846  1.96   2.005]]

Calculate the prediction

In [11]:
mus=priormean(xPred)
ypred=mus+np.dot(np.dot(Ks,KBInv),(yB-mxB))
print "Prediction: ",ypred
Prediction:  [  0.061   0.071   0.096   0.133   0.183   0.25    0.335   0.444   0.581
   0.75    0.951   1.182   1.439   1.715   2.003   2.299   2.599   2.905
   3.22    3.55    3.901   4.279   4.689   5.131   5.605   6.109   6.639
   7.191   7.764   8.354   8.961   9.583  10.223  10.88   11.554]

Calculate the covariance of the predictions

In [12]:
yvar=np.diag(Kss-np.dot(Ks,np.dot(KBInv,np.transpose(Ks))))
stds=np.sqrt(yvar)
print "Double Standard Deviation: ",2*stds
Double Standard Deviation:  [ 2.031  1.672  1.255  0.817  0.417  0.2    0.306  0.393  0.37   0.27   0.2
  0.257  0.328  0.328  0.257  0.2    0.27   0.37   0.393  0.306  0.2    0.417
  0.817  1.255  1.672  2.031  2.315  2.521  2.658  2.742  2.789  2.813
  2.824  2.829  2.831]

Plot training data and predicitons:

In [13]:
plt.figure(figsize=(12, 10))
plt.plot(x,mx,label="mean $m(x)$")
plt.hold(True)
plt.plot(xB,yB,'or',label="training data")
plt.plot(xPred,ypred,'--g',label="predictions")
plt.text(0.5,8,"$m(x)=0.25 \cdot x^2$ \n$k(x,x')=2 \cdot \exp(-0.5\cdot(x-x')^2)$ \n $\sigma_n^2=0.005$",fontsize=14)
plt.legend(loc=2,numpoints=1)
plt.title('Gaussian Process Prediction with prior quadratic mean')
plt.xlabel('x')
plt.ylabel('f(x)')
plt.axis([0,8,0,16])
# plot 2*standard deviation 95%-confidence interval
fillx = r_[xPred, xPred[::-1]]

filly = r_[ypred+2*stds, ypred[::-1]-2*stds[::-1]]
plt.fill(fillx, filly, facecolor='gray', edgecolor='white', alpha=0.3)
plt.show()