Merge branch 'master' of github.com:SheffieldML/GPy

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 <tuto_GP_regression.html>`_  
+* `A kernel overview <tuto_kernel_overview.html>`_ 
 * 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`
-

doc/tuto_GP_regression.rst

@@ -27,7 +27,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::
@@ -114,7 +114,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::

doc/tuto_kernel_overview.rst

@@ -0,0 +1,164 @@
+
+****************************
+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(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
+
+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. Here for example ``ker_add.D`` will return ``1`` whereas ``ker_prod.D`` will return ``2``.
+
+In order to add kernels defined on the different input spaces, the required command is::
+
+    ker_add_orth = k1.add_orthogonal(k2)
+
+.. figure::  Figures/tuto_kern_overview_add_orth.png
+    :align:  center
+    :height: 350px
+
+    Output of ``ker_add_orth.plot(plot_limits=[[-10,-10],[10,10]])``.
+
+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 part of the above equation. 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
+
+
+..  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')

grid_parameters.py

@@ -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()
-