.. _gp-tutorial: Gaussian Process ==================== Limbo relies on our C++-11 implementation of Gaussian processes (See :ref:`gaussian-process` for a short introduction ) which can be useful by itself. This tutorial explains how to create and use a Gaussian Process (GP). Data ---- We assume that our samples are in a vector called ``samples`` and that our observations are in a vector called ``observations``. Please note that the type of both observations and samples is Eigen::VectorXd (in this example, they are 1-D vectors). .. literalinclude:: ../../src/tutorials/gp.cpp :language: c++ :linenos: :lines: 75-85 Basic usage ------------ We first create a basic GP with an Exponential kernel (``kernel::Exp``) and a mean function equals to the mean of the obsevations (``mean::Data``). The ``Exp`` kernel needs a few parameters to be defined in a ``Params`` structure: .. literalinclude:: ../../src/tutorials/gp.cpp :language: c++ :linenos: :lines: 60-63 The type of the GP is defined by the following lines: .. literalinclude:: ../../src/tutorials/gp.cpp :language: c++ :linenos: :lines: 87-90 To use the GP, we need : - to instantiante a ``GP_t`` object - to call the method ``compute()``. .. literalinclude:: ../../src/tutorials/gp.cpp :language: c++ :linenos: :lines: 92-99 Here we assume that the noise is the same for all samples and that it is equal to 0.01. Querying the GP can be achieved in two different ways: - ``gp.mu(v)`` and ``gp.sigma(v)``, which return the mean and the variance (sigma squared) for the input data point ``v`` - ``std::tie(mu, sigma) = gp.query(v)``, which returns the mean and the variance at the same time. The second approach is faster because some computations are the same for ``mu`` and ``sigma``. To visualize the predictions of the GP, we can query it for many points and record the predictions in a file: .. literalinclude:: ../../src/tutorials/gp.cpp :language: c++ :linenos: :lines: 101-112 Hyper-parameter optimization ---------------------------- Most kernel functions have some parameters. It is common in the GP literature to set them by maximizing the log-likelihood of the data knowing the model (see :ref:`gaussian-process` for a description of this concept). In limbo, only a subset of the kernel functions can have their hyper-parameters optimized. The most common one is ``SquaredExpARD`` (Squared Exponential with Automatic Relevance Determination). A new GP type is defined as follows: .. literalinclude:: ../../src/tutorials/gp.cpp :language: c++ :linenos: :lines: 116-118 It uses the default values for the parameters of ``SquaredExpARD``: .. literalinclude:: ../../src/tutorials/gp.cpp :language: c++ :linenos: :lines: 64-65 After calling the ``compute()`` method, the hyper-parameters can be optimized by calling the ``optimize_hyperparams()`` function. The GP does not need to be recomputed and we pass ``false`` for the last parameter in ``compute()`` as we do not need to compute the kernel matrix again (it will be recomputed in the hyper-parameters optimization). .. literalinclude:: ../../src/tutorials/gp.cpp :language: c++ :linenos: :lines: 122-123 We can have a look at the difference between the two GPs: .. figure:: ../pics/gp.png :alt: Comparisons of GP :target: ../_images/gp.png This plot is generated using matplotlib: .. literalinclude:: ../../src/tutorials/plot_gp.py :language: python :linenos: