.. only:: html

    .. note::
        :class: sphx-glr-download-link-note

        Click :ref:`here <sphx_glr_download_auto_probabilistic_modeling_stochastic_processes_plot_gaussian_processes_comparison.py>`     to download the full example code
    .. rst-class:: sphx-glr-example-title

    .. _sphx_glr_auto_probabilistic_modeling_stochastic_processes_plot_gaussian_processes_comparison.py:


Comparison of covariance models for gaussian processes
======================================================

The main goal of this example is to briefly review the most important covariance models and compare them in terms of regularity. 

We first show how to define a covariance model, a temporal grid and a gaussian process. We first consider the squared exponential covariance model and show how the trajectories are sensitive to its parameters. We show how to define a trend. In the final section, we compare the trajectories from exponential and Matérn covariance models.

References
----------

* Carl Edward Rasmussen and Christopher K. I. Williams (2006) Gaussian Processes for Machine Learning. Chapter 4: "Covariance Functions", www.GaussianProcess.org/gpml

The anisotropic squared exponential model
-----------------------------------------

The `SquaredExponential` class allows to define covariance models :

* :math:`\sigma\in\mathbb{R}` is the amplitude parameter,
* :math:`\boldsymbol{\theta}\in\mathbb{R}^d` is the scale.


.. code-block:: default

    import openturns as ot
    import openturns.viewer as viewer
    from matplotlib import pylab as plt
    ot.Log.Show(ot.Log.NONE)








Amplitude values


.. code-block:: default

    amplitude = [3.5]
    # Scale values
    scale = [1.5]
    # Covariance model
    myModel = ot.SquaredExponential(scale, amplitude)








Gaussian processes
------------------

In order to create a `GaussianProcess`, we must have 
* a covariance model, 
* a grid. 

Optionnally, we can define a trend (we will see that later in the example). By default, the trend is zero. 

We consider the domain :math:`\mathcal{D}=[0,10]`. We discretize this domain with 100 cells (which corresponds to 101 nodes), with steps equal to 0.1 starting from 0: 

.. math::
   (x_0=x_{min}=0,\:x_1=0.1,\:\ldots,\:x_n=x_{max}=10).



.. code-block:: default

    xmin = 0.0
    step = 0.1
    n = 100
    myTimeGrid = ot.RegularGrid(xmin, step, n+1)
    graph = myTimeGrid.draw()
    graph.setTitle("Regular grid")
    view = viewer.View(graph)




.. image:: /auto_probabilistic_modeling/stochastic_processes/images/sphx_glr_plot_gaussian_processes_comparison_001.png
    :alt: Regular grid
    :class: sphx-glr-single-img





Then we create the gaussian process (by default the trend is zero). 


.. code-block:: default

    process = ot.GaussianProcess(myModel, myTimeGrid)
    process






.. raw:: html

    <p>GaussianProcess(trend=[x0]->[0.0], covariance=SquaredExponential(scale=[1.5], amplitude=[3.5]))</p>
    <br />
    <br />

Then we generate 10 trajectores with the `getSample` method. This trajectories are in a `ProcessSample`. 


.. code-block:: default

    nbTrajectories = 10
    sample = process.getSample(nbTrajectories)
    type(sample)








We can draw the trajectories with `drawMarginal`.


.. code-block:: default

    graph = sample.drawMarginal(0)
    graph.setTitle("amplitude=%.3f, scale=%.3f" % (amplitude[0], scale[0]))
    view = viewer.View(graph)





.. image:: /auto_probabilistic_modeling/stochastic_processes/images/sphx_glr_plot_gaussian_processes_comparison_002.png
    :alt: amplitude=3.500, scale=1.500
    :class: sphx-glr-single-img





In order to make the next examples easier, we define a function which plots a given number of trajectories from a gaussian process based on a covariance model. 


.. code-block:: default

    def plotCovarianceModel(myCovarianceModel,myTimeGrid,nbTrajectories):
        '''Plots the given number of trajectories with given covariance model.'''
        process = ot.GaussianProcess(myCovarianceModel, myTimeGrid)
        sample = process.getSample(nbTrajectories)
        graph = sample.drawMarginal(0)
        graph.setTitle("")
        return graph









The amplitude parameter sets the variance of the process. A greater amplitude increases the chances of getting larger absolute values of the process. 


.. code-block:: default

    amplitude = [7.]
    scale = [1.5]
    myModel = ot.SquaredExponential(scale, amplitude)
    graph = plotCovarianceModel(myModel,myTimeGrid,10)
    graph.setTitle("amplitude=%.3f, scale=%.3f" % (amplitude[0], scale[0]))
    view = viewer.View(graph)




.. image:: /auto_probabilistic_modeling/stochastic_processes/images/sphx_glr_plot_gaussian_processes_comparison_003.png
    :alt: amplitude=7.000, scale=1.500
    :class: sphx-glr-single-img





Modifying the scale parameter is here equivalent to stretch or contract the "time" :math:`x`. 


.. code-block:: default

    amplitude = [3.5]
    scale = [0.5]
    myModel = ot.SquaredExponential(scale, amplitude)
    graph = plotCovarianceModel(myModel,myTimeGrid,10)
    graph.setTitle("amplitude=%.3f, scale=%.3f" % (amplitude[0], scale[0]))
    view = viewer.View(graph)




.. image:: /auto_probabilistic_modeling/stochastic_processes/images/sphx_glr_plot_gaussian_processes_comparison_004.png
    :alt: amplitude=3.500, scale=0.500
    :class: sphx-glr-single-img





Define the trend
----------------

The trend is a deterministic function. With the `GaussianProcess` class, the associated process is the sum of a trend and a gaussian process with zero mean.


.. code-block:: default

    f = ot.SymbolicFunction(['x'], ['2*x'])
    fTrend = ot.TrendTransform(f,myTimeGrid)









.. code-block:: default

    amplitude = [3.5]
    scale = [1.5]
    myModel = ot.SquaredExponential(scale, amplitude)
    process = ot.GaussianProcess(fTrend,myModel, myTimeGrid)
    process






.. raw:: html

    <p>GaussianProcess(trend=[x]->[2*x], covariance=SquaredExponential(scale=[1.5], amplitude=[3.5]))</p>
    <br />
    <br />


.. code-block:: default

    nbTrajectories = 10
    sample = process.getSample(nbTrajectories)
    graph = sample.drawMarginal(0)
    graph.setTitle("amplitude=%.3f, scale=%.3f" % (amplitude[0], scale[0]))
    view = viewer.View(graph)




.. image:: /auto_probabilistic_modeling/stochastic_processes/images/sphx_glr_plot_gaussian_processes_comparison_005.png
    :alt: amplitude=3.500, scale=1.500
    :class: sphx-glr-single-img





Other covariance models
-----------------------

There are other covariance models. The models which are used more often are the following. 
* `SquaredExponential`. The generated processes can be derivated in mean square at all orders. 
* `MaternModel`. When :math:`\nu\rightarrow+\infty`, it converges to the squared exponential model. This model can be derivated :math:`k` times only if :math:`k<\nu`. In other words, when :math:`\nu` increases, then the trajectories are more and more regular. The particular case :math:`\nu=1/2` is the exponential model. The most commonly used values are :math:`\nu=3/2` and :math:`\nu=5/2`, which produce trajectories that are, in terms of regularity, in between the squared exponential and the exponential models. 
* `ExponentialModel`. The associated process is continus, but not differentiable.

The Matérn and exponential models
---------------------------------


.. code-block:: default

    amplitude = [1.0]
    scale = [1.0]
    nu1, nu2, nu3 = 2.5, 1.5, 0.5
    myModel1 = ot.MaternModel(scale, amplitude, nu1)
    myModel2 = ot.MaternModel(scale, amplitude, nu2)
    myModel3 = ot.MaternModel(scale, amplitude, nu3)









.. code-block:: default

    nbTrajectories = 10
    graph1 = plotCovarianceModel(myModel1,myTimeGrid,nbTrajectories)
    graph2 = plotCovarianceModel(myModel2,myTimeGrid,nbTrajectories)
    graph3 = plotCovarianceModel(myModel3,myTimeGrid,nbTrajectories)









.. code-block:: default

    from openturns.viewer import View
    import pylab as pl
    fig = pl.figure(figsize=(20, 6))
    ax1 = fig.add_subplot(1, 3, 1)
    _ = View(graph1, figure=fig, axes=[ax1])
    _ = ax1.set_title("Matern 5/2")
    ax2 = fig.add_subplot(1, 3, 2)
    _ = View(graph2, figure=fig, axes=[ax2])
    _ = ax2.set_title("Matern 3/2")
    ax3 = fig.add_subplot(1, 3, 3)
    _ = View(graph3, figure=fig, axes=[ax3])
    _ = ax3.set_title("Matern 1/2")




.. image:: /auto_probabilistic_modeling/stochastic_processes/images/sphx_glr_plot_gaussian_processes_comparison_006.png
    :alt: , Matern 5/2, Matern 3/2, Matern 1/2
    :class: sphx-glr-single-img





We see than, when :math:`\nu` increases, then the trajectories are smoother and smoother.


.. code-block:: default

    myExpModel = ot.ExponentialModel(scale, amplitude)









.. code-block:: default

    graph = plotCovarianceModel(myExpModel,myTimeGrid,nbTrajectories)
    graph.setTitle("Exponential")
    view = viewer.View(graph)




.. image:: /auto_probabilistic_modeling/stochastic_processes/images/sphx_glr_plot_gaussian_processes_comparison_007.png
    :alt: Exponential
    :class: sphx-glr-single-img





We see that the exponential model produces very irregular trajectories.


.. rst-class:: sphx-glr-timing

   **Total running time of the script:** ( 0 minutes  0.815 seconds)


.. _sphx_glr_download_auto_probabilistic_modeling_stochastic_processes_plot_gaussian_processes_comparison.py:


.. only :: html

 .. container:: sphx-glr-footer
    :class: sphx-glr-footer-example



  .. container:: sphx-glr-download sphx-glr-download-python

     :download:`Download Python source code: plot_gaussian_processes_comparison.py <plot_gaussian_processes_comparison.py>`



  .. container:: sphx-glr-download sphx-glr-download-jupyter

     :download:`Download Jupyter notebook: plot_gaussian_processes_comparison.ipynb <plot_gaussian_processes_comparison.ipynb>`


.. only:: html

 .. rst-class:: sphx-glr-signature

    `Gallery generated by Sphinx-Gallery <https://sphinx-gallery.github.io>`_