diff --git a/doc/Figures/kern-def.png b/doc/Figures/kern-def.png new file mode 100644 index 00000000..bad43b09 Binary files /dev/null and b/doc/Figures/kern-def.png differ diff --git a/doc/Figures/tuto_GP_regression_m1.png b/doc/Figures/tuto_GP_regression_m1.png index e4174825..e1a11fb1 100644 Binary files a/doc/Figures/tuto_GP_regression_m1.png and b/doc/Figures/tuto_GP_regression_m1.png differ diff --git a/doc/Figures/tuto_GP_regression_m2.png b/doc/Figures/tuto_GP_regression_m2.png index bf28f4e0..7e54e919 100644 Binary files a/doc/Figures/tuto_GP_regression_m2.png and b/doc/Figures/tuto_GP_regression_m2.png differ diff --git a/doc/Figures/tuto_kern_overview_add_orth.png b/doc/Figures/tuto_kern_overview_add_orth.png new file mode 100644 index 00000000..0d4f1c4e Binary files /dev/null and b/doc/Figures/tuto_kern_overview_add_orth.png differ diff --git a/doc/Figures/tuto_kern_overview_allkern.png b/doc/Figures/tuto_kern_overview_allkern.png new file mode 100644 index 00000000..f3406b07 Binary files /dev/null and b/doc/Figures/tuto_kern_overview_allkern.png differ diff --git a/doc/Figures/tuto_kern_overview_basicdef.png b/doc/Figures/tuto_kern_overview_basicdef.png new file mode 100644 index 00000000..bad43b09 Binary files /dev/null and b/doc/Figures/tuto_kern_overview_basicdef.png differ diff --git a/doc/Figures/tuto_kern_overview_mANOVA.png b/doc/Figures/tuto_kern_overview_mANOVA.png new file mode 100644 index 00000000..db49e3bd Binary files /dev/null and b/doc/Figures/tuto_kern_overview_mANOVA.png differ diff --git a/doc/Figures/tuto_kern_overview_mANOVAdec.png b/doc/Figures/tuto_kern_overview_mANOVAdec.png new file mode 100644 index 00000000..ef154263 Binary files /dev/null and b/doc/Figures/tuto_kern_overview_mANOVAdec.png differ diff --git a/doc/index.rst b/doc/index.rst index 28690e99..b62ff6a7 100644 --- a/doc/index.rst +++ b/doc/index.rst @@ -8,8 +8,8 @@ Welcome to GPy's documentation! For a quick start, you can have a look at one of the tutorials: * `Basic Gaussian process regression `_ +* `A kernel overview `_ * Advanced GP regression (Forthcoming) -* Kernel manipulation (Forthcoming) * Writting kernels (Forthcoming) You may also be interested by some examples in the GPy/examples folder. @@ -28,4 +28,3 @@ Indices and tables * :ref:`genindex` * :ref:`modindex` * :ref:`search` - diff --git a/doc/tuto_GP_regression.rst b/doc/tuto_GP_regression.rst index 3527e86f..f14afe62 100644 --- a/doc/tuto_GP_regression.rst +++ b/doc/tuto_GP_regression.rst @@ -12,7 +12,7 @@ We first import the libraries we will need: :: import numpy as np import GPy -1 dimensional model +1-dimensional model =================== For this toy example, we assume we have the following inputs and outputs:: @@ -99,7 +99,7 @@ Once again, we can use ``print(m)`` and ``m.plot()`` to look at the resulting mo GP regression model after optimization of the parameters. -2 dimensional example +2-dimensional example ===================== Here is a 2 dimensional example:: diff --git a/doc/tuto_kernel_overview.rst b/doc/tuto_kernel_overview.rst new file mode 100644 index 00000000..9ef86589 --- /dev/null +++ b/doc/tuto_kernel_overview.rst @@ -0,0 +1,169 @@ + +**************************** +tutorial : A kernel overview +**************************** + +First we import the libraries we will need :: + + import pylab as pb + import numpy as np + import GPy + pb.ion() + +For most kernels, the dimension is the only mandatory parameter to define a kernel object. However, it is also possible to specify the values of the parameters. For example, the three following commands are valid for defining a squared exponential kernel (ie rbf or Gaussian) :: + + ker1 = GPy.kern.rbf(D=1) # Equivalent to ker1 = GPy.kern.rbf(D=1, variance=1., lengthscale=1.) + ker2 = GPy.kern.rbf(D=1, variance = 1.5, lengthscale=2.) + ker3 = GPy.kern.rbf(1, .5, .5) + +A `plot` and a `print` functions are implemented to represent kernel objects :: + + print ker1 + + ker1.plot() + ker2.plot() + ker3.plot() + +.. figure:: Figures/tuto_kern_overview_basicdef.png + :align: center + :height: 350px + +:: + + import pylab as pb + import numpy as np + import GPy + pb.ion() + + ker1 = GPy.kern.rbf(D=1) # Equivalent to ker1 = GPy.kern.rbf(D=1, variance=1., lengthscale=1.) + ker2 = GPy.kern.rbf(D=1, variance = .75, lengthscale=3.) + ker3 = GPy.kern.rbf(1, .5, .25) + + ker1.plot() + ker2.plot() + ker3.plot() + #pb.savefig("Figures/tuto_kern_overview_basicdef.png") + + kernels = [GPy.kern.rbf(1), GPy.kern.exponential(1), GPy.kern.Matern32(1), GPy.kern.Matern52(1), GPy.kern.Brownian(1), GPy.kern.bias(1), GPy.kern.linear(1), GPy.kern.spline(1), GPy.kern.periodic_exponential(1), GPy.kern.periodic_Matern32(1), GPy.kern.periodic_Matern52(1), GPy.kern.white(1)] + kernel_names = ["GPy.kern.rbf", "GPy.kern.exponential", "GPy.kern.Matern32", "GPy.kern.Matern52", "GPy.kern.Brownian", "GPy.kern.bias", "GPy.kern.linear", "GPy.kern.spline", "GPy.kern.periodic_exponential", "GPy.kern.periodic_Matern32", "GPy.kern.periodic_Matern52", "GPy.kern.white"] + + pb.figure(figsize=(16,12)) + pb.subplots_adjust(wspace=.5, hspace=.5) + for i, kern in enumerate(kernels): + pb.subplot(3,4,i+1) + kern.plot(x=7.5,plot_limits=[0.00001,15.]) + pb.title(kernel_names[i]+ '\n') + #pb.axes([.1,.1,.8,.7]) + #pb.figtext(.5,.9,'Foo Bar', fontsize=18, ha='center') + #pb.figtext(.5,.85,'Lorem ipsum dolor sit amet, consectetur adipiscing elit',fontsize=10,ha='center') + + # actual plot for the noise + i = 11 + X = np.linspace(0.,15.,201) + WN = 0*X + WN[100] = 1. + pb.subplot(3,4,i+1) + pb.plot(X,WN,'b') + +Implemented kernels +=================== + +Many kernels are already implemented in GPy. Here is a summary of most of them: + +.. figure:: Figures/tuto_kern_overview_allkern.png + :align: center + :height: 800px + +On the other hand, it is possible to use the `sympy` package to build new kernels. This will be the subject of another tutorial. + +Operations to combine kernel +============================ + +In ``GPy``, kernel objects can be combined with the usual ``+`` and ``*`` operators. :: + + k1 = GPy.kern.rbf(1,variance=1., lengthscale=2) + k2 = GPy.kern.Matern32(1,variance=1., lengthscale=2) + + ker_add = k1 + k2 + print ker_add + + ker_prod = k1 * k2 + print ker_prod + +Note that by default, the operator ``+`` adds kernels defined on the same input space whereas ``*`` assumes that the kernels are defined on different input spaces. :: + + ker_add.D + ker_prod.D + +In order to add kernels defined on the different input spaces, the required command is:: + + ker_add_orth = k1.add_orthogonal(k2) + +The resulting kernel is + ker_add_orth.plot(plot_limits=[[-10,-10],[10,10]]) + +.. figure:: Figures/tuto_kern_overview_add_orth.png + :align: center + :height: 350px + +Example : Building an ANOVA kernel +================================== + +In two dimensions ANOVA kernels have the following form: + +.. math:: + + k_{ANOVA}(x,y) = \prod_{i=1}^2 (1 + k_i(x_i,y_i)) = 1 + k_1(x_1,y_1) + k_2(x_2,y_2) + k_1(x_1,y_1) \times k_2(x_2,y_2). + +Let us assume that we want to define an ANOVA kernel with a Matern 3/2 kernel for :math:`k_i`. As seen previously, we can define this kernel as follow:: + + k_cst = GPy.kern.bias(1,variance=1.) + k_mat = GPy.kern.Matern52(1,variance=1., lengthscale=3) + Kanova = (k_cst + k_mat) * (k_cst + k_mat) + print Kanova + +Note the ties between the lengthscales of ``Kanova`` to keep the number of lengthscales equal to 2. On the other hand, there are four variance terms in the new parameterization: one for each term of the right hand sign of the equation above. We can illustrate the use of this kernel on a toy example:: + + # sample inputs and outputs + X = np.random.uniform(-3.,3.,(40,2)) + Y = 0.5*X[:,:1] + 0.5*X[:,1:] + 2*np.sin(X[:,:1]) * np.sin(X[:,1:]) + + # Create GP regression model + m = GPy.models.GP_regression(X,Y,Kanova) + m.plot() + + +.. figure:: Figures/tuto_kern_overview_mANOVA.png + :align: center + :height: 350px + +As :math:`k_{ANOVA}` corresponds to the sum of 4 kernels, the best predictor can be splited in a sum of 4 functions + +.. math:: + + bp(x) & = k(x)^t K^{-1} Y \\ + & = (1 + k_1(x_1) + k_2(x_2) + k_1(x_1)k_2(x_2))^t K^{-1} Y \\ + & = 1^t K^{-1} Y + k_1(x_1)^t K^{-1} Y + k_2(x_2)^t K^{-1} Y + (k_1(x_1)k_2(x_2))^t K^{-1} Y + +The submodels can be represented with the option ``which_function`` of ``plot``: :: + + pb.figure(figsize=(20,5)) + pb.subplots_adjust(wspace=0.5) + pb.subplot(1,5,1) + m.plot() + pb.subplot(1,5,2) + pb.ylabel("= ",rotation='horizontal',fontsize='30') + pb.subplot(1,5,3) + m.plot(which_functions=[False,True,False,False]) + pb.ylabel("cst +",rotation='horizontal',fontsize='30') + pb.subplot(1,5,4) + m.plot(which_functions=[False,False,True,False]) + pb.ylabel("+ ",rotation='horizontal',fontsize='30') + pb.subplot(1,5,5) + pb.ylabel("+ ",rotation='horizontal',fontsize='30') + m.plot(which_functions=[False,False,False,True]) + + +.. figure:: Figures/tuto_kern_overview_mANOVAdec.png + :align: center + :height: 200px diff --git a/grid_parameters.py b/grid_parameters.py deleted file mode 100644 index 64d82755..00000000 --- a/grid_parameters.py +++ /dev/null @@ -1,64 +0,0 @@ -import numpy as np -import pylab as pb -pb.ion() -import sys -import GPy - -pb.close('all') - -N = 200 -M = 15 -resolution=5 - -X = np.linspace(0,12,N)[:,None] -Z = np.linspace(0,12,M)[:,None] # inducing points (fixed for now) -Y = np.sin(X) + np.random.randn(*X.shape)/np.sqrt(50.) -#k = GPy.kern.rbf(1) -k = GPy.kern.Matern32(1) + GPy.kern.white(1) - -models = [GPy.models.sparse_GP_regression(X,Y,Z=Z,kernel=k) - ,GPy.models.sparse_GP_regression(X,Y,Z=Z,kernel=k) - ,GPy.models.sparse_GP_regression(X,Y,Z=Z,kernel=k) - ,GPy.models.sparse_GP_regression(X,Y,Z=Z,kernel=k)] -models[0].scale_factor = 1. -models[1].scale_factor = 10. -models[2].scale_factor = 100. -models[3].scale_factor = 1000. - #GPy.models.sgp_debugB(X,Y,Z=Z,kernel=k), - #GPy.models.sgp_debugC(X,Y,Z=Z,kernel=k)]#, - #GPy.models.sgp_debugE(X,Y,Z=Z,kernel=k)] - -[m.constrain_fixed('white',0.1) for m in models] - -#xx,yy = np.mgrid[1.5:4:0+resolution*1j,-2:2:0+resolution*1j] -xx,yy = np.mgrid[3:16:0+resolution*1j,-2:1:0+resolution*1j] - -lls = [] -cgs = [] -grads = [] -count = 0 -for l,v in zip(xx.flatten(),yy.flatten()): - count += 1 - print count, 'of', resolution**2 - sys.stdout.flush() - - [m.set('lengthscale',l) for m in models] - [m.set('_variance',10.**v) for m in models] - lls.append([m.log_likelihood() for m in models]) - grads.append([m.log_likelihood_gradients() for m in models]) - cgs.append([m.checkgrad(verbose=0,return_ratio=True) for m in models]) - -lls = np.array(zip(*lls)).reshape(-1,resolution,resolution) -cgs = np.array(zip(*cgs)).reshape(-1,resolution,resolution) - -for ll,cg in zip(lls,cgs): - pb.figure() - pb.contourf(xx,yy,ll,100,cmap=pb.cm.gray) - pb.colorbar() - try: - pb.contour(xx,yy,np.exp(ll),colors='k') - except: - pass - pb.scatter(xx.flatten(),yy.flatten(),20,np.log(np.abs(cg.flatten())),cmap=pb.cm.jet,linewidth=0) - pb.colorbar() -