<SYSTEM_TASK:>
Given the following text description, write Python code to implement the functionality described below step by step
<END_TASK>
<USER_TASK:>
Description:
 Lab session 1
Step1:  1 Covariance Functions
GPy supports a wide range of covariance functions. Each function is suitable for some applications. In this section, we work with the covariance functions in GPy. Let's start with defining an squared exponential or rbf covariance function in one dimension
Step2:  A summary of the kernel can be obtained using the command print k.
Step3:  It is also possible to plot the kernel as a function of one of its inputs (whilst fixing the other) with k.plot(). Use "?" to view the properties of the plot function for kern instance.
Step4:  Setting Covariance Function Parameters
The value of the covariance function parameters can be accessed and modified using k.* where the * can refer to the parameter name as it appears in print k. We'll now set the lengthscale of the covariance to different values, and then plot the resulting covariance using the k.plot() method.
Step5:  Exercise 1
a) What is the effect of the lengthscale parameter on the covariance function?
b) Now change the code used above for plotting the covariances associated with the length scale to see the influence of the variance parameter. What is the effect of the the variance parameter on the covariance function?
Step6:  c) Instead of rbf, try constructing and plotting the following  covariance functions
Step7:  Computing the kernel matrix given input data, $\mathbf{X}$
Let $\mathbf{X}$ be a $n$ × $d$ numpy array. Given a kernel $k$, the covariance matrix associated to
$\mathbf{X}$ is obtained with C = k.K(X,X) . The positive semi-definiteness of $k$ ensures that C
is a positive semi-definite (psd) matrix regardless of the initial points $\mathbf{X}$. This can be
checked numerically by looking at the eigenvalues
Step8:  Combining Covariance Functions
In GPy you can easily combine covariance functions you have created using the sum and product operators, + and *. For example, we can combine an exponentiated quadratic covariance with a Matern 5/2 as follows
Step9:  It is also possible to multiply two kernel functions
Step10:  2 Sampling from a Gaussian Process
Gaussian process provides a prior over an infinite dimensional function. It is defined by a covariance function and a mean function. When we compute the covariance matrix using kern.K(X, X) i.e. compute a covariance matrix between the values of the function that correspond to the input locations in the matrix X. Using this we can draw sample paths from a Gaussian process
Step11:  We can see the structure of the covariance matrix we are plotting from if we visualize C.
Step12:  Exercise 2
a) Try a range of different covariance functions and values to plot the corresponding sample paths for each using the same approach given above.
Step13:  b) Can you tell the covariance structures that have been used for generating the
sample paths shown in the figure below?
<br>
<center>
<img src="http
Step14:  A GP regression model is defined by first specifying the covariance function for analysis. Then an instance of the model is generated with a default set of parameters. Then it is possible to view the parameters using print m and visualize the posterior mean prediction and variances using m.plot
Step15:  The actual predictions of the model for a set of points Xstar can be computed using m.predict(Xstar)
Step16:  Exercise 3
a) What do you think about this first fit? Does the prior given by the GP seem to be
adapted?
b) The parameters of the models can be modified using a regular expression matching the parameters names (for example m['noise'] = 0.001 ). Change the values of the parameters to obtain a better fit.
Step17:  c) Random sample paths from the conditional GP can be obtained using np.random.multivariate_normal(mu[
Step18:  Covariance Function Parameter Estimation
The parameters values can be estimated by maximizing the likelihood of the observations. Since we don’t want one of the variance to become negative during the optimization, we can constrain all parameters to be positive before running the optimisation.
Step19:  We can optimize the hyperparameters of the model using the m.optimize() method.
Step20:  Exercise 4
a) Generate random function samples using the optmized model and plot alongside the input data
Step21:  b) Modify the kernel used for building the model to investigate its influence on the model
Step22:  4 A Running Example
Now we will consider a small example with real world data
Step23:  Exercise 5
a) Build a Gaussian process model for the olympic data set using a combination of an RBF and a bias covariance function. Fit the covariance function parameters and the noise to the data. Plot the fit and error bars from 1870 to 2030. Do you think the predictions are reasonable? If not why not?
Step24:  b) Fit the same model, but this time intialize the length scale of the RBF kernel to 0.5. What has happened? Which of model has the higher log likelihood, this one or the one from (a)? 
Hint
Step25:  c) Modify your model by including two covariance functions. Intitialize a covariance function with an exponentiated quadratic part, a Matern 3/2 part and a bias covariance. Set the initial lengthscale of the exponentiated quadratic to 80 years, set the initial length scale of the Matern 3/2 to 10 years. Optimize the new model and plot the fit again. How does it compare with the previous model?
Step26:  d) Repeat part c) but now initialize both of the covariance functions' lengthscales to 20 years. Check the model parameters, what happens now?
Step27:  e) Now model the data with a product of an exponentiated quadratic covariance function and a linear covariance function. Fit the covariance function parameters. Why are the variance parameters of the linear part so small? How could this be fixed?
<ASSISTANT_TASK:>
Python Code:
# import python modules
import GPy
import numpy as np
from matplotlib import pyplot as plt
# call matplotlib with the inline command to make plots appear within the browser
%matplotlib inline
Explanation: Lab session 1: Gaussian Process models with GPy
Source: Gaussian Process Summer School 2015
The aim of this lab session is to get us started with GPy Library. The current draft of the online documentation of GPy is available from this page. We will focus on three aspects of GPs: the kernel, the random sample paths and the GP regression model.
Requirements:
* GPy: Installation instructions available on the homepage.
* Scipy Stack: This includes numpy, matplotlib and Ipython. Installation can be done using pip:
(sudo) pip install numpy --upgrade
(sudo) pip install jupyter --upgrade
(sudo) pip install matplotlib --upgrade
(sudo) pip install ipython[all] --upgrade
The sudo is optional if you want to have installation in the root folder when working in the Linux operating system. It should not be used for Anaconda or when working in the Windows operating system.
Anaconda: Necessary for Windows, optional for Linux.
For Windows Operating System:
* Install Anaconda by downloading from this link.
* Install GPy by opening the Command Prompt window and typing the following command:
pip install GPy
* All the dependencies required for running the code are either available in Anaconda or installed using GPy.
End of explanation
# The documentation to use the RBF function. There are several advanced options such as useGPU which are 
# important for practical applications. The "?" symbol can be used with any function or class to view its
# documentation
GPy.kern.RBF?
# input dimension
d = 1          
# variance
var = 1.       
# lengthscale
length = 0.2
# define the kernel
k = GPy.kern.RBF(d, variance=var, lengthscale=length)
Explanation: 1 Covariance Functions
GPy supports a wide range of covariance functions. Each function is suitable for some applications. In this section, we work with the covariance functions in GPy. Let's start with defining an squared exponential or rbf covariance function in one dimension:
$$
k(\mathbf{x},\mathbf{x'}) = \sigma^2 \text{exp} \left( - \frac{\|\mathbf{x} - \mathbf{x'} \|}{2l^2} \right)
$$
End of explanation
# view the parameters of the covariance function
print k
Explanation: A summary of the kernel can be obtained using the command print k.
End of explanation
# plot the covariance function
k.plot()
Explanation: It is also possible to plot the kernel as a function of one of its inputs (whilst fixing the other) with k.plot(). Use "?" to view the properties of the plot function for kern instance.
End of explanation
# by default, all the parameters are set to 1. for the RBF kernel
k = GPy.kern.RBF(d)     
# we experiment with different length scale parameter values here
theta = np.asarray([0.2,0.5,1.,2.,4.])
# create an instance of a figure
fig = plt.figure()
ax = plt.subplot(111)
# iterate over the lengthscales
for t in theta:
    k.lengthscale=t
    # plot in the same figure with a different color
    k.plot(ax=ax, color=np.random.rand(3,), plot_limits=[-10.0,10.0])
plt.legend(theta)
Explanation: Setting Covariance Function Parameters
The value of the covariance function parameters can be accessed and modified using k.* where the * can refer to the parameter name as it appears in print k. We'll now set the lengthscale of the covariance to different values, and then plot the resulting covariance using the k.plot() method.
End of explanation
# by default, all the parameters are set to 1. for the RBF kernel
k = GPy.kern.RBF(d)     
# we experiment with different length scale parameter values here
var = np.asarray([0.2,0.5,1.,2.,4.])
# create an instance of a figure
fig = plt.figure()
ax = plt.subplot(111)
# iterate over the lengthscales
for v in var:
    k.variance = v
    # plot in the same figure with a different color
    k.plot(ax=ax, color=np.random.rand(3,))
plt.legend(var)
Explanation: Exercise 1
a) What is the effect of the lengthscale parameter on the covariance function?
b) Now change the code used above for plotting the covariances associated with the length scale to see the influence of the variance parameter. What is the effect of the the variance parameter on the covariance function?
End of explanation
# look for the kernel documentation
GPy.kern.Matern32?
# input dim
d = 1
# create the Matern32 kernel
k = GPy.kern.Matern32(d)
# view the kernel
print k
k.plot()
Explanation: c) Instead of rbf, try constructing and plotting the following  covariance functions: exponential, Matern32, Matern52, Brownian, linear, bias, rbfcos, periodic_Matern32, etc. Use the tab key to look for the functions GPy.kern submodule or search in the GPy documentation.
End of explanation
# input data: 50*2 matrix of iid standard Gaussians
X = np.random.rand(50,2)       
# create the matern52 kernel
k = GPy.kern.Matern52(input_dim=2)
# compute the kernel matrix
C = k.K(X,X)
# computes eigenvalues of matrix
eigvals = np.linalg.eigvals(C) 
# plot the eigen values
plt.bar(np.arange(len(eigvals)), eigvals)
plt.title('Eigenvalues of the Matern 5/2 Covariance')
Explanation: Computing the kernel matrix given input data, $\mathbf{X}$
Let $\mathbf{X}$ be a $n$ × $d$ numpy array. Given a kernel $k$, the covariance matrix associated to
$\mathbf{X}$ is obtained with C = k.K(X,X) . The positive semi-definiteness of $k$ ensures that C
is a positive semi-definite (psd) matrix regardless of the initial points $\mathbf{X}$. This can be
checked numerically by looking at the eigenvalues:
End of explanation
# define rbf and matern52 kernels
kern1 = GPy.kern.RBF(1, variance=1., lengthscale=2.)
kern2 = GPy.kern.Matern52(1, variance=2., lengthscale=4.)
# combine both kernels
kern = kern1 + kern2
print kern
kern.plot(plot_limits=[-7,9])
Explanation: Combining Covariance Functions
In GPy you can easily combine covariance functions you have created using the sum and product operators, + and *. For example, we can combine an exponentiated quadratic covariance with a Matern 5/2 as follows:
End of explanation
kern = kern1*kern2
print kern
kern.plot(plot_limits=[-6,8])
Explanation: It is also possible to multiply two kernel functions:
End of explanation
# define RBF kernel
k = GPy.kern.RBF(input_dim=1,lengthscale=0.2)
# define X to be 500 points evenly spaced over [0,1]
X = np.linspace(0.,1.,500) 
# make the numpy array to 2D array 
X = X[:,None] 
# set mean function i.e. 0 everywhere
mu = np.zeros((500)) 
# compute covariance matrix associated with inputs X
C = k.K(X,X) 
# Generate 20 separate samples paths from a Gaussian with mean mu and covariance C
Z = np.random.multivariate_normal(mu,C,20)
# open a new plotting window
fig = plt.figure()     
for i in range(20):
    plt.plot(X[:],Z[i,:])
Explanation: 2 Sampling from a Gaussian Process
Gaussian process provides a prior over an infinite dimensional function. It is defined by a covariance function and a mean function. When we compute the covariance matrix using kern.K(X, X) i.e. compute a covariance matrix between the values of the function that correspond to the input locations in the matrix X. Using this we can draw sample paths from a Gaussian process:
$$
\mathbf{f} = \mathcal{N}(\mathbf{X}|\mathbf{0},C)
$$
End of explanation
plt.matshow(C)
Explanation: We can see the structure of the covariance matrix we are plotting from if we visualize C.
End of explanation
# Define input points and mean function
X = np.atleast_2d(np.linspace(0.,1.,500)).T
mu = np.zeros((500)) 
# sample paths for RBF kernel with different parameters
k = GPy.kern.RBF(input_dim=1,lengthscale=0.05)
C = k.K(X,X) 
# Generate 20 separate samples paths from a Gaussian with mean mu and covariance C
Z = np.random.multivariate_normal(mu,C,10)
# open a new plotting window
fig = plt.figure()     
for i in range(10):
    plt.plot(X[:],Z[i,:])
# sample paths for Matern32 kernel with different parameters
k = GPy.kern.PeriodicMatern52(input_dim=1, lengthscale=0.01, period=1)
C = k.K(X,X) 
# Generate 20 separate samples paths from a Gaussian with mean mu and covariance C
Z = np.random.multivariate_normal(mu,C,10)
# open a new plotting window
fig = plt.figure()     
for i in range(10):
    plt.plot(X[:],Z[i,:])
Explanation: Exercise 2
a) Try a range of different covariance functions and values to plot the corresponding sample paths for each using the same approach given above.
End of explanation
# input data points
X = np.linspace(0.05,0.95,10)[:,None]
# generate observations through function f
Y = -np.cos(np.pi*X) + np.sin(4*np.pi*X) + np.random.normal(loc=0.0, scale=0.1, size=(10,1)) 
# plot the generated data
plt.figure()
plt.plot(X,Y,'kx',mew=1.5)
Explanation: b) Can you tell the covariance structures that have been used for generating the
sample paths shown in the figure below?
<br>
<center>
<img src="http://ml.dcs.shef.ac.uk/gpss/gpws14/figa.png" alt="Figure a" style="width: 30%;"> 
<img src="http://ml.dcs.shef.ac.uk/gpss/gpws14/figb.png" alt="Figure b" style="width: 30%;"> 
<img src="http://ml.dcs.shef.ac.uk/gpss/gpws14/figd.png" alt="Figure d" style="width: 30%;">
</center>
3 Gaussian Process Regression Model
We will combine the Gaussian process prior with data to form a GP regression model with GPy. We will generate data from the function $f ( x ) = − \cos(\pi x ) + \sin(4\pi x )$ over $[0, 1]$, adding some noise to give $y(x) = f(x) + \epsilon$, with the noise being Gaussian distributed, $\epsilon \sim \mathcal{N}(0, 0.01)$.
End of explanation
# create instance of kernel
k = GPy.kern.RBF(input_dim=1, variance=1., lengthscale=1.)
# create instance of GP regression model
m = GPy.models.GPRegression(X,Y,k)
# view model parameters
print m
# visualize posterior mean and variances
m.plot()
Explanation: A GP regression model is defined by first specifying the covariance function for analysis. Then an instance of the model is generated with a default set of parameters. Then it is possible to view the parameters using print m and visualize the posterior mean prediction and variances using m.plot:
End of explanation
# obtain 5 test points
Xstar = np.linspace(0.01,0.99,5)[:,None]
# predict the output for the test points
Ystar, Vstar = m.predict(Xstar)
# print results
print Ystar
Explanation: The actual predictions of the model for a set of points Xstar can be computed using m.predict(Xstar):
End of explanation
# set the noise parameter of the model using desired SNR
SNR = 10.0
m.Gaussian_noise.variance = m.rbf.variance/SNR
# check the model parameters and plot the model
print m
m.plot()
Explanation: Exercise 3
a) What do you think about this first fit? Does the prior given by the GP seem to be
adapted?
b) The parameters of the models can be modified using a regular expression matching the parameters names (for example m['noise'] = 0.001 ). Change the values of the parameters to obtain a better fit.
End of explanation
# get desired input points where model needs to be evaluated
Xp = np.linspace(0.0,1.0,100)[:,None]
# obtain posterior mean and variances
mu, C = m.predict(Xp, full_cov=True)
# generate 10 random paths of the distribution
nPaths = 10
paths = np.random.multivariate_normal(mu[:,0], C, nPaths)
# plot the dataset and the generated paths
plt.figure()
plt.plot(X,Y,'kx',mew=5)
for i in range(nPaths):
    plt.plot(Xp[:],paths[i,:])
Explanation: c) Random sample paths from the conditional GP can be obtained using np.random.multivariate_normal(mu[:,0],C) where the mean vector and covariance matrix mu, C are obtained through the predict function mu, C = m.predict(Xp,full_cov=True). Obtain 10 samples from the posterior sample and plot them alongside the data below.
End of explanation
m.constrain_positive()
Explanation: Covariance Function Parameter Estimation
The parameters values can be estimated by maximizing the likelihood of the observations. Since we don’t want one of the variance to become negative during the optimization, we can constrain all parameters to be positive before running the optimisation.
End of explanation
m.optimize()
m.plot()
print m
Explanation: We can optimize the hyperparameters of the model using the m.optimize() method.
End of explanation
# get desired input points where model needs to be evaluated
Xp = np.linspace(0.0,1.0,100)[:,None]
# obtain posterior mean and variances
mu, C = m.predict(Xp, full_cov=True)
# generate 10 random paths of the distribution
nPaths = 10
paths = np.random.multivariate_normal(mu[:,0], C, nPaths)
# plot the dataset and the generated paths
plt.figure()
plt.plot(X,Y,'kx',mew=5)
for i in range(nPaths):
    plt.plot(Xp[:],paths[i,:])
Explanation: Exercise 4
a) Generate random function samples using the optmized model and plot alongside the input data
End of explanation
# create instance of kernel and GP Regression model
k = GPy.kern.Matern32(input_dim=1, variance=1., lengthscale=1.)
m = GPy.models.GPRegression(X,Y,k)
# optimize the model
m.constrain_positive()
m.optimize()
# view model parameters
print m
m.plot()
# get desired input points where model needs to be evaluated
Xp = np.linspace(0.0,1.0,100)[:,None]
# obtain posterior mean and variances
mu, C = m.predict(Xp, full_cov=True)
# generate 10 random paths of the distribution
nPaths = 10
paths = np.random.multivariate_normal(mu[:,0], C, nPaths)
# plot the dataset and the generated paths
plt.figure()
plt.plot(X,Y,'kx',mew=5)
for i in range(nPaths):
    plt.plot(Xp[:],paths[i,:])
Explanation: b) Modify the kernel used for building the model to investigate its influence on the model
End of explanation
# load marathon timing dataset
data = np.genfromtxt('marathon.csv', delimiter=',')
# set the input and output data
X = data[:, 0:1]
Y = data[:, 1:2]
# plot the timings
plt.plot(X, Y, 'bx')
plt.xlabel('year')
plt.ylabel('marathon pace min/km')
Explanation: 4 A Running Example
Now we will consider a small example with real world data: data giving the pace of all marathons run at the olympics:
End of explanation
# create the covariance function
kern = GPy.kern.RBF(1) + GPy.kern.Bias(1)
# create the GP Regression model
model = GPy.models.GPRegression(X, Y, kern)
# optimize the model
model.optimize()
model.plot()
Explanation: Exercise 5
a) Build a Gaussian process model for the olympic data set using a combination of an RBF and a bias covariance function. Fit the covariance function parameters and the noise to the data. Plot the fit and error bars from 1870 to 2030. Do you think the predictions are reasonable? If not why not?
End of explanation
# create the covariance function and set the desired parameter values
kern2 = GPy.kern.RBF(1, variance=0.5) + GPy.kern.Bias(1)
model2 = GPy.models.GPRegression(X, Y, kern2)
# optimize and plot the model
model2.optimize()
model2.plot()
# compare the log likelihoods
print model.log_likelihood(), model2.log_likelihood()
Explanation: b) Fit the same model, but this time intialize the length scale of the RBF kernel to 0.5. What has happened? Which of model has the higher log likelihood, this one or the one from (a)? 
Hint: use model.log_likelihood() for computing the log likelihood.
End of explanation
# create the covariance function and set the desired parameter values
kern3 = GPy.kern.RBF(1, lengthscale=80.0) + GPy.kern.Matern32(1, lengthscale=10.0) + GPy.kern.Bias(1)
model3 = GPy.models.GPRegression(X, Y, kern3)
# optimize and plot the model
model3.optimize()
model3.plot()
# compare the log likelihoods
print model.log_likelihood(), model3.log_likelihood()
Explanation: c) Modify your model by including two covariance functions. Intitialize a covariance function with an exponentiated quadratic part, a Matern 3/2 part and a bias covariance. Set the initial lengthscale of the exponentiated quadratic to 80 years, set the initial length scale of the Matern 3/2 to 10 years. Optimize the new model and plot the fit again. How does it compare with the previous model?
End of explanation
# create the covariance function and set the desired parameter values
kern4 = GPy.kern.RBF(1, lengthscale=20.0) + GPy.kern.Matern32(1, lengthscale=20.0) + GPy.kern.Bias(1)
model4 = GPy.models.GPRegression(X, Y, kern4)
# optimize and plot the model
model4.optimize()
model4.plot()
# compare the log likelihoods
print model3.log_likelihood(), model4.log_likelihood()
Explanation: d) Repeat part c) but now initialize both of the covariance functions' lengthscales to 20 years. Check the model parameters, what happens now?
End of explanation
# create the covariance function and set the desired parameter values
kern5 = GPy.kern.RBF(1, lengthscale=5.0, variance=5.0)*GPy.kern.Linear(1)
model5 = GPy.models.GPRegression(X, Y, kern5)
# optimize and plot the model
model5.optimize()
model5.plot()
# get the model parameters
print model5
# compare the log likelihoods
print model.log_likelihood(), model5.log_likelihood()
Explanation: e) Now model the data with a product of an exponentiated quadratic covariance function and a linear covariance function. Fit the covariance function parameters. Why are the variance parameters of the linear part so small? How could this be fixed?
End of explanation
<END_TASK>