{ "cells": [ { "cell_type": "code", "execution_count": 1, "metadata": { "collapsed": true }, "outputs": [], "source": [ "import GPy" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "As of Mon 12th of Oct running on devel branch of GPy 0.8.8" ] }, { "cell_type": "code", "execution_count": 2, "metadata": { "collapsed": false }, "outputs": [], "source": [ "GPy.plotting.change_plotting_library('plotly')" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "# Gaussian process regression tutorial\n", "\n", "### Nicolas Durrande 2013\n", "#### with edits by James Hensman and Neil D. Lawrence\n", "\n", "We will see in this tutorial the basics for building a 1 dimensional and a 2 dimensional Gaussian process regression model, also known as a kriging model.\n", "\n", "We first import the libraries we will need:" ] }, { "cell_type": "code", "execution_count": 3, "metadata": { "collapsed": false }, "outputs": [], "source": [ "import numpy as np" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "## 1-dimensional model\n", "\n", "For this toy example, we assume we have the following inputs and outputs:" ] }, { "cell_type": "code", "execution_count": 4, "metadata": { "collapsed": false }, "outputs": [], "source": [ "X = np.random.uniform(-3.,3.,(20,1))\n", "Y = np.sin(X) + np.random.randn(20,1)*0.05" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "Note that the observations Y include some noise.\n", "\n", "The first step is to define the covariance kernel we want to use for the model. We choose here a kernel based on Gaussian kernel (i.e. rbf or square exponential):" ] }, { "cell_type": "code", "execution_count": 5, "metadata": { "collapsed": false }, "outputs": [], "source": [ "kernel = GPy.kern.RBF(input_dim=1, variance=1., lengthscale=1.)" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "The parameter input_dim stands for the dimension of the input space. The parameters `variance` and `lengthscale` are optional, and default to 1. Many other kernels are implemented, type `GPy.kern.` to see a list" ] }, { "cell_type": "code", "execution_count": 6, "metadata": { "collapsed": false }, "outputs": [], "source": [ "#type GPy.kern. here:\n", "GPy.kern.BasisFuncKernel?" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "The inputs required for building the model are the observations and the kernel:" ] }, { "cell_type": "code", "execution_count": 7, "metadata": { "collapsed": false }, "outputs": [], "source": [ "m = GPy.models.GPRegression(X,Y,kernel)" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "By default, some observation noise is added to the model. The functions `display` and `plot` give an insight of the model we have just built:" ] }, { "cell_type": "code", "execution_count": 8, "metadata": { "collapsed": false }, "outputs": [ { "data": { "text/html": [ "\n", "\n", "

\n", "Model: GP regression
\n", "Objective: 22.9717924697
\n", "Number of Parameters: 3
\n", "Number of Optimization Parameters: 3
\n", "Updates: True
\n", "

\n", "\n", "\n", "\n", "\n", "\n", "
GP_regression. valueconstraintspriors
rbf.variance 1.0 +ve
rbf.lengthscale 1.0 +ve
Gaussian_noise.variance 1.0 +ve
" ], "text/plain": [ "" ] }, "metadata": {}, "output_type": "display_data" } ], "source": [ "from IPython.display import display\n", "display(m)" ] }, { "cell_type": "code", "execution_count": 9, "metadata": { "collapsed": false }, "outputs": [ { "name": "stdout", "output_type": "stream", "text": [ "This is the format of your plot grid:\n", "[ (1,1) x1,y1 ]\n", "\n" ] }, { "data": { "text/html": [ "" ], "text/plain": [ "" ] }, "execution_count": 9, "metadata": {}, "output_type": "execute_result" } ], "source": [ "fig = m.plot()\n", "GPy.plotting.show(fig, filename='basic_gp_regression_notebook')" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "The above cell shows our GP regression model before optimization of the parameters. The shaded region corresponds to ~95% confidence intervals (ie +/- 2 standard deviation).\n", "\n", "The default values of the kernel parameters may not be optimal for the current data (for example, the confidence intervals seems too wide on the previous figure). A common approach is to find the values of the parameters that maximize the likelihood of the data. It as easy as calling `m.optimize` in GPy:" ] }, { "cell_type": "code", "execution_count": 10, "metadata": { "collapsed": false }, "outputs": [], "source": [ "m.optimize(messages=True)" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "If we want to perform some restarts to try to improve the result of the optimization, we can use the `optimize_restarts` function. This selects random (drawn from $N(0,1)$) initializations for the parameter values, optimizes each, and sets the model to the best solution found." ] }, { "cell_type": "code", "execution_count": 11, "metadata": { "collapsed": false }, "outputs": [ { "name": "stdout", "output_type": "stream", "text": [ "Optimization restart 1/10, f = -15.1436482683\n", "Optimization restart 2/10, f = -15.1436482683\n", "Optimization restart 3/10, f = -15.1436482682\n", "Optimization restart 4/10, f = -15.1436482682\n", "Optimization restart 5/10, f = -15.1436482682\n", "Optimization restart 6/10, f = -15.1436482682\n", "Optimization restart 7/10, f = -15.1436482683\n", "Optimization restart 8/10, f = -15.1436482682\n", "Optimization restart 9/10, f = -15.1436482683\n", "Optimization restart 10/10, f = -15.1436482683\n" ] }, { "data": { "text/plain": [ "[,\n", " ,\n", " ,\n", " ,\n", " ,\n", " ,\n", " ,\n", " ,\n", " ,\n", " ,\n", " ]" ] }, "execution_count": 11, "metadata": {}, "output_type": "execute_result" } ], "source": [ "m.optimize_restarts(num_restarts = 10)" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "In this simple example, the objective function (usually!) has only one local minima, and each of the found solutions are the same. \n", "\n", "Once again, we can use `print(m)` and `m.plot()` to look at the resulting model resulting model. This time, the paraemters values have been optimized agains the log likelihood (aka the log marginal likelihood): the fit shoul dbe much better. " ] }, { "cell_type": "code", "execution_count": 12, "metadata": { "collapsed": false }, "outputs": [ { "data": { "text/html": [ "\n", "\n", "

\n", "Model: GP regression
\n", "Objective: -15.1436482683
\n", "Number of Parameters: 3
\n", "Number of Optimization Parameters: 3
\n", "Updates: True
\n", "

\n", "\n", "\n", "\n", "\n", "\n", "
GP_regression. valueconstraintspriors
rbf.variance 1.35354271667 +ve
rbf.lengthscale 1.94630136743 +ve
Gaussian_noise.variance0.00248112830273 +ve
" ], "text/plain": [ "" ] }, "metadata": {}, "output_type": "display_data" }, { "name": "stdout", "output_type": "stream", "text": [ "This is the format of your plot grid:\n", "[ (1,1) x1,y1 ]\n", "\n" ] }, { "data": { "text/html": [ "" ], "text/plain": [ "" ] }, "execution_count": 12, "metadata": {}, "output_type": "execute_result" } ], "source": [ "display(m)\n", "fig = m.plot()\n", "GPy.plotting.show(fig, filename='basic_gp_regression_notebook_optimized')" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "### New plotting of GPy 0.9 and later\n", "The new plotting allows you to plot the density of a GP object more fine grained by plotting more percentiles of the distribution color coded by their opacity" ] }, { "cell_type": "code", "execution_count": 13, "metadata": { "collapsed": false }, "outputs": [ { "data": { "text/html": [ "\n", "\n", "

\n", "Model: GP regression
\n", "Objective: -15.1436482683
\n", "Number of Parameters: 3
\n", "Number of Optimization Parameters: 3
\n", "Updates: True
\n", "

\n", "\n", "\n", "\n", "\n", "\n", "
GP_regression. valueconstraintspriors
rbf.variance 1.35354271667 +ve
rbf.lengthscale 1.94630136743 +ve
Gaussian_noise.variance0.00248112830273 +ve
" ], "text/plain": [ "" ] }, "metadata": {}, "output_type": "display_data" }, { "name": "stdout", "output_type": "stream", "text": [ "This is the format of your plot grid:\n", "[ (1,1) x1,y1 ]\n", "\n" ] }, { "data": { "text/html": [ "" ], "text/plain": [ "" ] }, "execution_count": 13, "metadata": {}, "output_type": "execute_result" } ], "source": [ "display(m)\n", "fig = m.plot(plot_density=True)\n", "GPy.plotting.show(fig, filename='basic_gp_regression_density_notebook_optimized')" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "## 2-dimensional example\n", "\n", "Here is a 2 dimensional example:" ] }, { "cell_type": "code", "execution_count": 14, "metadata": { "collapsed": false }, "outputs": [ { "name": "stdout", "output_type": "stream", "text": [ "This is the format of your plot grid:\n", "[ (1,1) x1,y1 ]\n", "\n" ] }, { "data": { "text/html": [ "" ], "text/plain": [ "" ] }, "metadata": {}, "output_type": "display_data" }, { "data": { "text/html": [ "\n", "\n", "

\n", "Model: GP regression
\n", "Objective: -24.7900663215
\n", "Number of Parameters: 5
\n", "Number of Optimization Parameters: 5
\n", "Updates: True
\n", "

\n", "\n", "\n", "\n", "\n", "\n", "\n", "
GP_regression. valueconstraintspriors
sum.Mat52.variance 0.361421808902 +ve
sum.Mat52.lengthscale (2,) +ve
sum.white.variance 0.000644606566433 +ve
Gaussian_noise.variance0.000644606566433 +ve
" ], "text/plain": [ "" ] }, "metadata": {}, "output_type": "display_data" } ], "source": [ "# sample inputs and outputs\n", "X = np.random.uniform(-3.,3.,(50,2))\n", "Y = np.sin(X[:,0:1]) * np.sin(X[:,1:2])+np.random.randn(50,1)*0.05\n", "\n", "# define kernel\n", "ker = GPy.kern.Matern52(2,ARD=True) + GPy.kern.White(2)\n", "\n", "# create simple GP model\n", "m = GPy.models.GPRegression(X,Y,ker)\n", "\n", "# optimize and plot\n", "m.optimize(messages=True,max_f_eval = 1000)\n", "fig = m.plot()\n", "display(GPy.plotting.show(fig, filename='basic_gp_regression_notebook_2d'))\n", "display(m)" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "The flag `ARD=True` in the definition of the `Matern` kernel specifies that we want one lengthscale parameter per dimension (ie the GP is not isotropic). Note that for 2-d plotting, only the mean is shown. " ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "## Plotting slices\n", "To see the uncertaintly associated with the above predictions, we can plot slices through the surface. this is done by passing the optional `fixed_inputs` argument to the plot function. `fixed_inputs` is a list of tuples containing which of the inputs to fix, and to which value.\n", "\n", "To get horixontal slices of the above GP, we'll fix second (index 1) input to -1, 0, and 1.5:" ] }, { "cell_type": "code", "execution_count": 15, "metadata": { "collapsed": false }, "outputs": [ { "name": "stdout", "output_type": "stream", "text": [ "This is the format of your plot grid:\n", "[ (1,1) x1,y1 ]\n", "[ (2,1) x1,y2 ]\n", "[ (3,1) x1,y3 ]\n", "\n" ] }, { "data": { "text/html": [ "" ], "text/plain": [ "" ] }, "execution_count": 15, "metadata": {}, "output_type": "execute_result" } ], "source": [ "slices = [-1, 0, 1.5]\n", "figure = GPy.plotting.plotting_library().figure(3, 1, \n", " shared_xaxes=True,\n", " subplot_titles=('slice at -1', \n", " 'slice at 0', \n", " 'slice at 1.5', \n", " )\n", " )\n", "for i, y in zip(range(3), slices):\n", " canvas = m.plot(figure=figure, fixed_inputs=[(1,y)], row=(i+1), plot_data=False)\n", "GPy.plotting.show(canvas, filename='basic_gp_regression_notebook_slicing')" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "A few things to note:\n", " * we've also passed the optional `ax` argument, to mnake the GP plot on a particular subplot\n", " * the data look strange here: we're seeing slices of the GP, but all the data are displayed, even though they might not be close to the current slice." ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "To get vertical slices, we simply fixed the other input. We'll turn the display of data off also:" ] }, { "cell_type": "code", "execution_count": 16, "metadata": { "collapsed": false }, "outputs": [ { "name": "stdout", "output_type": "stream", "text": [ "This is the format of your plot grid:\n", "[ (1,1) x1,y1 ]\n", "[ (2,1) x1,y2 ]\n", "[ (3,1) x1,y3 ]\n", "\n" ] }, { "data": { "text/html": [ "" ], "text/plain": [ "" ] }, "execution_count": 16, "metadata": {}, "output_type": "execute_result" } ], "source": [ "slices = [-1, 0, 1.5]\n", "figure = GPy.plotting.plotting_library().figure(3, 1, \n", " shared_xaxes=True,\n", " subplot_titles=('slice at -1', \n", " 'slice at 0', \n", " 'slice at 1.5', \n", " )\n", " )\n", "for i, y in zip(range(3), slices):\n", " canvas = m.plot(figure=figure, fixed_inputs=[(0,y)], row=(i+1), plot_data=False)\n", "GPy.plotting.show(canvas, filename='basic_gp_regression_notebook_slicing_vertical')" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "You can find a host of other plotting options in the `m.plot` docstring. `Type m.plot?` to see. " ] } ], "metadata": { "kernelspec": { "display_name": "Python 2", "language": "python", "name": "python2" }, "language_info": { "codemirror_mode": { "name": "ipython", "version": 2 }, "file_extension": ".py", "mimetype": "text/x-python", "name": "python", "nbconvert_exporter": "python", "pygments_lexer": "ipython2", "version": "2.7.11" } }, "nbformat": 4, "nbformat_minor": 0 }