.. DO NOT EDIT. .. THIS FILE WAS AUTOMATICALLY GENERATED BY SPHINX-GALLERY. .. TO MAKE CHANGES, EDIT THE SOURCE PYTHON FILE: .. "auto_probabilistic_modeling/stochastic_processes/plot_gaussian_processes_comparison.py" .. LINE NUMBERS ARE GIVEN BELOW. .. only:: html .. note:: :class: sphx-glr-download-link-note :ref:`Go to the end ` to download the full example code .. rst-class:: sphx-glr-example-title .. _sphx_glr_auto_probabilistic_modeling_stochastic_processes_plot_gaussian_processes_comparison.py: Compare covariance models ========================= .. GENERATED FROM PYTHON SOURCE LINES 6-11 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. .. GENERATED FROM PYTHON SOURCE LINES 13-17 References ---------- * Carl Edward Rasmussen and Christopher K. I. Williams (2006) Gaussian Processes for Machine Learning. Chapter 4: "Covariance Functions", www.GaussianProcess.org/gpml .. GENERATED FROM PYTHON SOURCE LINES 19-26 The anisotropic squared exponential model ----------------------------------------- The :class:`~openturns.SquaredExponential` class allows one to define covariance models: * :math:`\sigma\in\mathbb{R}` is the amplitude parameter, * :math:`\boldsymbol{\theta}\in\mathbb{R}^d` is the scale. .. GENERATED FROM PYTHON SOURCE LINES 28-35 .. code-block:: Python import pylab as pl from openturns.viewer import View import openturns as ot import openturns.viewer as viewer ot.Log.Show(ot.Log.NONE) .. GENERATED FROM PYTHON SOURCE LINES 36-37 Amplitude values .. GENERATED FROM PYTHON SOURCE LINES 37-43 .. code-block:: Python amplitude = [3.5] # Scale values scale = [1.5] # Covariance model myModel = ot.SquaredExponential(scale, amplitude) .. GENERATED FROM PYTHON SOURCE LINES 44-59 Gaussian processes ------------------ In order to create a :class:`~openturns.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). .. GENERATED FROM PYTHON SOURCE LINES 61-69 .. code-block:: Python 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-sg:: /auto_probabilistic_modeling/stochastic_processes/images/sphx_glr_plot_gaussian_processes_comparison_001.png :alt: Regular grid :srcset: /auto_probabilistic_modeling/stochastic_processes/images/sphx_glr_plot_gaussian_processes_comparison_001.png :class: sphx-glr-single-img .. GENERATED FROM PYTHON SOURCE LINES 70-71 Then we create the gaussian process (by default the trend is zero). .. GENERATED FROM PYTHON SOURCE LINES 73-76 .. code-block:: Python process = ot.GaussianProcess(myModel, myTimeGrid) process .. raw:: html
class=GaussianProcess mesh=class=Mesh name=Unnamed dimension=1 vertices=class=Sample name=Unnamed implementation=class=SampleImplementation name=Unnamed size=101 dimension=1 description=[t] data=[[0],[0.1],[0.2],[0.3],[0.4],[0.5],[0.6],[0.7],[0.8],[0.9],[1],[1.1],[1.2],[1.3],[1.4],[1.5],[1.6],[1.7],[1.8],[1.9],[2],[2.1],[2.2],[2.3],[2.4],[2.5],[2.6],[2.7],[2.8],[2.9],[3],[3.1],[3.2],[3.3],[3.4],[3.5],[3.6],[3.7],[3.8],[3.9],[4],[4.1],[4.2],[4.3],[4.4],[4.5],[4.6],[4.7],[4.8],[4.9],[5],[5.1],[5.2],[5.3],[5.4],[5.5],[5.6],[5.7],[5.8],[5.9],[6],[6.1],[6.2],[6.3],[6.4],[6.5],[6.6],[6.7],[6.8],[6.9],[7],[7.1],[7.2],[7.3],[7.4],[7.5],[7.6],[7.7],[7.8],[7.9],[8],[8.1],[8.2],[8.3],[8.4],[8.5],[8.6],[8.7],[8.8],[8.9],[9],[9.1],[9.2],[9.3],[9.4],[9.5],[9.6],[9.7],[9.8],[9.9],[10]] simplices=[[0,1],[1,2],[2,3],[3,4],[4,5],[5,6],[6,7],[7,8],[8,9],[9,10],[10,11],[11,12],[12,13],[13,14],[14,15],[15,16],[16,17],[17,18],[18,19],[19,20],[20,21],[21,22],[22,23],[23,24],[24,25],[25,26],[26,27],[27,28],[28,29],[29,30],[30,31],[31,32],[32,33],[33,34],[34,35],[35,36],[36,37],[37,38],[38,39],[39,40],[40,41],[41,42],[42,43],[43,44],[44,45],[45,46],[46,47],[47,48],[48,49],[49,50],[50,51],[51,52],[52,53],[53,54],[54,55],[55,56],[56,57],[57,58],[58,59],[59,60],[60,61],[61,62],[62,63],[63,64],[64,65],[65,66],[66,67],[67,68],[68,69],[69,70],[70,71],[71,72],[72,73],[73,74],[74,75],[75,76],[76,77],[77,78],[78,79],[79,80],[80,81],[81,82],[82,83],[83,84],[84,85],[85,86],[86,87],[87,88],[88,89],[89,90],[90,91],[91,92],[92,93],[93,94],[94,95],[95,96],[96,97],[97,98],[98,99],[99,100]] trend=class=TrendTransform inherited from class=VertexValueFunction evaluation=class=Function name=Unnamed implementation=class=FunctionImplementation name=Unnamed description=[x0,x0,y0] evaluationImplementation=class=TrendEvaluation name=Unnamed function=class=Function name=Unnamed implementation=class=FunctionImplementation name=Unnamed description=[x0,y0] evaluationImplementation=class=SymbolicEvaluation name=Unnamed inputVariablesNames=[x0] outputVariablesNames=[y0] formulas=[0.0] gradientImplementation=class=SymbolicGradient name=Unnamed evaluation=class=SymbolicEvaluation name=Unnamed inputVariablesNames=[x0] outputVariablesNames=[y0] formulas=[0.0] hessianImplementation=class=SymbolicHessian name=Unnamed evaluation=class=SymbolicEvaluation name=Unnamed inputVariablesNames=[x0] outputVariablesNames=[y0] formulas=[0.0] gradientImplementation=class=CenteredFiniteDifferenceGradient name=Unnamed epsilon=class=Point name=Unnamed dimension=2 values=[1e-05,1e-05] evaluation=class=TrendEvaluation name=Unnamed function=class=Function name=Unnamed implementation=class=FunctionImplementation name=Unnamed description=[x0,y0] evaluationImplementation=class=SymbolicEvaluation name=Unnamed inputVariablesNames=[x0] outputVariablesNames=[y0] formulas=[0.0] gradientImplementation=class=SymbolicGradient name=Unnamed evaluation=class=SymbolicEvaluation name=Unnamed inputVariablesNames=[x0] outputVariablesNames=[y0] formulas=[0.0] hessianImplementation=class=SymbolicHessian name=Unnamed evaluation=class=SymbolicEvaluation name=Unnamed inputVariablesNames=[x0] outputVariablesNames=[y0] formulas=[0.0] hessianImplementation=class=CenteredFiniteDifferenceHessian name=Unnamed epsilon=class=Point name=Unnamed dimension=2 values=[0.0001,0.0001] evaluation=class=TrendEvaluation name=Unnamed function=class=Function name=Unnamed implementation=class=FunctionImplementation name=Unnamed description=[x0,y0] evaluationImplementation=class=SymbolicEvaluation name=Unnamed inputVariablesNames=[x0] outputVariablesNames=[y0] formulas=[0.0] gradientImplementation=class=SymbolicGradient name=Unnamed evaluation=class=SymbolicEvaluation name=Unnamed inputVariablesNames=[x0] outputVariablesNames=[y0] formulas=[0.0] hessianImplementation=class=SymbolicHessian name=Unnamed evaluation=class=SymbolicEvaluation name=Unnamed inputVariablesNames=[x0] outputVariablesNames=[y0] formulas=[0.0] covarianceModel=class=SquaredExponential scale=class=Point name=Unnamed dimension=1 values=[1.5] amplitude=class=Point name=Unnamed dimension=1 values=[3.5] covarianceCholeskyFactor=class=TriangularMatrix dimension=0 implementation=class=MatrixImplementation name=Unnamed rows=0 columns=0 values=[] isInitialized=false hasStationaryTrend=true checkedStationaryTrend=true


.. GENERATED FROM PYTHON SOURCE LINES 77-78 Then we generate 10 trajectores with the `getSample` method. This trajectories are in a `ProcessSample`. .. GENERATED FROM PYTHON SOURCE LINES 80-84 .. code-block:: Python nbTrajectories = 10 sample = process.getSample(nbTrajectories) type(sample) .. GENERATED FROM PYTHON SOURCE LINES 85-86 We can draw the trajectories with `drawMarginal`. .. GENERATED FROM PYTHON SOURCE LINES 88-93 .. code-block:: Python graph = sample.drawMarginal(0) graph.setTitle("amplitude=%.3f, scale=%.3f" % (amplitude[0], scale[0])) view = viewer.View(graph) .. image-sg:: /auto_probabilistic_modeling/stochastic_processes/images/sphx_glr_plot_gaussian_processes_comparison_002.png :alt: amplitude=3.500, scale=1.500 :srcset: /auto_probabilistic_modeling/stochastic_processes/images/sphx_glr_plot_gaussian_processes_comparison_002.png :class: sphx-glr-single-img .. GENERATED FROM PYTHON SOURCE LINES 94-95 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. .. GENERATED FROM PYTHON SOURCE LINES 98-107 .. code-block:: Python 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 .. GENERATED FROM PYTHON SOURCE LINES 108-109 The amplitude parameter sets the variance of the process. A greater amplitude increases the chances of getting larger absolute values of the process. .. GENERATED FROM PYTHON SOURCE LINES 111-118 .. code-block:: Python amplitude = [7.0] 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-sg:: /auto_probabilistic_modeling/stochastic_processes/images/sphx_glr_plot_gaussian_processes_comparison_003.png :alt: amplitude=7.000, scale=1.500 :srcset: /auto_probabilistic_modeling/stochastic_processes/images/sphx_glr_plot_gaussian_processes_comparison_003.png :class: sphx-glr-single-img .. GENERATED FROM PYTHON SOURCE LINES 119-120 Modifying the scale parameter is here equivalent to stretch or contract the "time" :math:`x`. .. GENERATED FROM PYTHON SOURCE LINES 122-129 .. code-block:: Python 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-sg:: /auto_probabilistic_modeling/stochastic_processes/images/sphx_glr_plot_gaussian_processes_comparison_004.png :alt: amplitude=3.500, scale=0.500 :srcset: /auto_probabilistic_modeling/stochastic_processes/images/sphx_glr_plot_gaussian_processes_comparison_004.png :class: sphx-glr-single-img .. GENERATED FROM PYTHON SOURCE LINES 130-134 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. .. GENERATED FROM PYTHON SOURCE LINES 136-139 .. code-block:: Python f = ot.SymbolicFunction(["x"], ["2*x"]) fTrend = ot.TrendTransform(f, myTimeGrid) .. GENERATED FROM PYTHON SOURCE LINES 140-146 .. code-block:: Python amplitude = [3.5] scale = [1.5] myModel = ot.SquaredExponential(scale, amplitude) process = ot.GaussianProcess(fTrend, myModel, myTimeGrid) process .. raw:: html
class=GaussianProcess mesh=class=Mesh name=Unnamed dimension=1 vertices=class=Sample name=Unnamed implementation=class=SampleImplementation name=Unnamed size=101 dimension=1 description=[t] data=[[0],[0.1],[0.2],[0.3],[0.4],[0.5],[0.6],[0.7],[0.8],[0.9],[1],[1.1],[1.2],[1.3],[1.4],[1.5],[1.6],[1.7],[1.8],[1.9],[2],[2.1],[2.2],[2.3],[2.4],[2.5],[2.6],[2.7],[2.8],[2.9],[3],[3.1],[3.2],[3.3],[3.4],[3.5],[3.6],[3.7],[3.8],[3.9],[4],[4.1],[4.2],[4.3],[4.4],[4.5],[4.6],[4.7],[4.8],[4.9],[5],[5.1],[5.2],[5.3],[5.4],[5.5],[5.6],[5.7],[5.8],[5.9],[6],[6.1],[6.2],[6.3],[6.4],[6.5],[6.6],[6.7],[6.8],[6.9],[7],[7.1],[7.2],[7.3],[7.4],[7.5],[7.6],[7.7],[7.8],[7.9],[8],[8.1],[8.2],[8.3],[8.4],[8.5],[8.6],[8.7],[8.8],[8.9],[9],[9.1],[9.2],[9.3],[9.4],[9.5],[9.6],[9.7],[9.8],[9.9],[10]] simplices=[[0,1],[1,2],[2,3],[3,4],[4,5],[5,6],[6,7],[7,8],[8,9],[9,10],[10,11],[11,12],[12,13],[13,14],[14,15],[15,16],[16,17],[17,18],[18,19],[19,20],[20,21],[21,22],[22,23],[23,24],[24,25],[25,26],[26,27],[27,28],[28,29],[29,30],[30,31],[31,32],[32,33],[33,34],[34,35],[35,36],[36,37],[37,38],[38,39],[39,40],[40,41],[41,42],[42,43],[43,44],[44,45],[45,46],[46,47],[47,48],[48,49],[49,50],[50,51],[51,52],[52,53],[53,54],[54,55],[55,56],[56,57],[57,58],[58,59],[59,60],[60,61],[61,62],[62,63],[63,64],[64,65],[65,66],[66,67],[67,68],[68,69],[69,70],[70,71],[71,72],[72,73],[73,74],[74,75],[75,76],[76,77],[77,78],[78,79],[79,80],[80,81],[81,82],[82,83],[83,84],[84,85],[85,86],[86,87],[87,88],[88,89],[89,90],[90,91],[91,92],[92,93],[93,94],[94,95],[95,96],[96,97],[97,98],[98,99],[99,100]] trend=class=TrendTransform inherited from class=VertexValueFunction evaluation=class=Function name=Unnamed implementation=class=FunctionImplementation name=Unnamed description=[x,x0,y0] evaluationImplementation=class=TrendEvaluation name=Unnamed function=class=Function name=Unnamed implementation=class=FunctionImplementation name=Unnamed description=[x,y0] evaluationImplementation=class=SymbolicEvaluation name=Unnamed inputVariablesNames=[x] outputVariablesNames=[y0] formulas=[2*x] gradientImplementation=class=SymbolicGradient name=Unnamed evaluation=class=SymbolicEvaluation name=Unnamed inputVariablesNames=[x] outputVariablesNames=[y0] formulas=[2*x] hessianImplementation=class=SymbolicHessian name=Unnamed evaluation=class=SymbolicEvaluation name=Unnamed inputVariablesNames=[x] outputVariablesNames=[y0] formulas=[2*x] gradientImplementation=class=CenteredFiniteDifferenceGradient name=Unnamed epsilon=class=Point name=Unnamed dimension=2 values=[1e-05,1e-05] evaluation=class=TrendEvaluation name=Unnamed function=class=Function name=Unnamed implementation=class=FunctionImplementation name=Unnamed description=[x,y0] evaluationImplementation=class=SymbolicEvaluation name=Unnamed inputVariablesNames=[x] outputVariablesNames=[y0] formulas=[2*x] gradientImplementation=class=SymbolicGradient name=Unnamed evaluation=class=SymbolicEvaluation name=Unnamed inputVariablesNames=[x] outputVariablesNames=[y0] formulas=[2*x] hessianImplementation=class=SymbolicHessian name=Unnamed evaluation=class=SymbolicEvaluation name=Unnamed inputVariablesNames=[x] outputVariablesNames=[y0] formulas=[2*x] hessianImplementation=class=CenteredFiniteDifferenceHessian name=Unnamed epsilon=class=Point name=Unnamed dimension=2 values=[0.0001,0.0001] evaluation=class=TrendEvaluation name=Unnamed function=class=Function name=Unnamed implementation=class=FunctionImplementation name=Unnamed description=[x,y0] evaluationImplementation=class=SymbolicEvaluation name=Unnamed inputVariablesNames=[x] outputVariablesNames=[y0] formulas=[2*x] gradientImplementation=class=SymbolicGradient name=Unnamed evaluation=class=SymbolicEvaluation name=Unnamed inputVariablesNames=[x] outputVariablesNames=[y0] formulas=[2*x] hessianImplementation=class=SymbolicHessian name=Unnamed evaluation=class=SymbolicEvaluation name=Unnamed inputVariablesNames=[x] outputVariablesNames=[y0] formulas=[2*x] covarianceModel=class=SquaredExponential scale=class=Point name=Unnamed dimension=1 values=[1.5] amplitude=class=Point name=Unnamed dimension=1 values=[3.5] covarianceCholeskyFactor=class=TriangularMatrix dimension=0 implementation=class=MatrixImplementation name=Unnamed rows=0 columns=0 values=[] isInitialized=false hasStationaryTrend=false checkedStationaryTrend=false


.. GENERATED FROM PYTHON SOURCE LINES 147-153 .. code-block:: Python 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-sg:: /auto_probabilistic_modeling/stochastic_processes/images/sphx_glr_plot_gaussian_processes_comparison_005.png :alt: amplitude=3.500, scale=1.500 :srcset: /auto_probabilistic_modeling/stochastic_processes/images/sphx_glr_plot_gaussian_processes_comparison_005.png :class: sphx-glr-single-img .. GENERATED FROM PYTHON SOURCE LINES 154-168 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 continuous, but not differentiable. .. GENERATED FROM PYTHON SOURCE LINES 170-172 The Matérn and exponential models --------------------------------- .. GENERATED FROM PYTHON SOURCE LINES 174-181 .. code-block:: Python 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) .. GENERATED FROM PYTHON SOURCE LINES 182-187 .. code-block:: Python nbTrajectories = 10 graph1 = plotCovarianceModel(myModel1, myTimeGrid, nbTrajectories) graph2 = plotCovarianceModel(myModel2, myTimeGrid, nbTrajectories) graph3 = plotCovarianceModel(myModel3, myTimeGrid, nbTrajectories) .. GENERATED FROM PYTHON SOURCE LINES 188-199 .. code-block:: Python 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-sg:: /auto_probabilistic_modeling/stochastic_processes/images/sphx_glr_plot_gaussian_processes_comparison_006.png :alt: , Matern 5/2, Matern 3/2, Matern 1/2 :srcset: /auto_probabilistic_modeling/stochastic_processes/images/sphx_glr_plot_gaussian_processes_comparison_006.png :class: sphx-glr-single-img .. GENERATED FROM PYTHON SOURCE LINES 200-201 We see than, when :math:`\nu` increases, then the trajectories are smoother and smoother. .. GENERATED FROM PYTHON SOURCE LINES 203-205 .. code-block:: Python myExpModel = ot.ExponentialModel(scale, amplitude) .. GENERATED FROM PYTHON SOURCE LINES 206-210 .. code-block:: Python graph = plotCovarianceModel(myExpModel, myTimeGrid, nbTrajectories) graph.setTitle("Exponential") view = viewer.View(graph) .. image-sg:: /auto_probabilistic_modeling/stochastic_processes/images/sphx_glr_plot_gaussian_processes_comparison_007.png :alt: Exponential :srcset: /auto_probabilistic_modeling/stochastic_processes/images/sphx_glr_plot_gaussian_processes_comparison_007.png :class: sphx-glr-single-img .. GENERATED FROM PYTHON SOURCE LINES 211-212 We see that the exponential model produces very irregular trajectories. .. _sphx_glr_download_auto_probabilistic_modeling_stochastic_processes_plot_gaussian_processes_comparison.py: .. only:: html .. container:: sphx-glr-footer sphx-glr-footer-example .. container:: sphx-glr-download sphx-glr-download-jupyter :download:`Download Jupyter notebook: plot_gaussian_processes_comparison.ipynb ` .. container:: sphx-glr-download sphx-glr-download-python :download:`Download Python source code: plot_gaussian_processes_comparison.py `