""" Compare covariance models ========================= """ # %% # The main goal of this example is to briefly review the most important covariance models and compare them in terms of regularity of the trajectories. # # 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 :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. # %% import openturns.viewer as otv import openturns as ot import matplotlib.pyplot as plt # %% # Amplitude values amplitude = [3.5] # Scale values scale = [1.5] # Covariance model myModel = ot.SquaredExponential(scale, amplitude) # %% # 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). # # %% xmin = 0.0 step = 0.1 n = 100 myTimeGrid = ot.RegularGrid(xmin, step, n + 1) graph = myTimeGrid.draw() graph.setTitle("Regular grid") view = otv.View(graph) # %% # Then we create the Gaussian process (by default the trend is zero). # %% process = ot.GaussianProcess(myModel, myTimeGrid) # %% # Then we generate 10 trajectores with the `getSample` method. This trajectories are in a :class:`~openturns.ProcessSample`. # %% nbTrajectories = 10 sample = process.getSample(nbTrajectories) type(sample) # %% # We can draw the trajectories with `drawMarginal`. # %% graph = sample.drawMarginal(0) graph.setTitle("amplitude=%.3f, scale=%.3f" % (amplitude[0], scale[0])) view = otv.View(graph) # %% # 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. # %% 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. # %% 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 = otv.View(graph) # %% # Modifying the scale parameter is here equivalent to stretch or contract the "time" :math:`x`. # %% 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 = otv.View(graph) # %% # Define the trend # ---------------- # # The trend is a deterministic function. With the :class:`~openturns.GaussianProcess` class, the associated process is the sum of a trend and a Gaussian process with zero mean. # %% f = ot.SymbolicFunction(["x"], ["2*x"]) fTrend = ot.TrendTransform(f, myTimeGrid) # %% amplitude = [3.5] scale = [1.5] myModel = ot.SquaredExponential(scale, amplitude) process = ot.GaussianProcess(fTrend, myModel, myTimeGrid) # %% # sphinx_gallery_thumbnail_number = 5 nbTrajectories = 10 sample = process.getSample(nbTrajectories) graph = sample.drawMarginal(0) graph.setTitle("amplitude=%.3f, scale=%.3f" % (amplitude[0], scale[0])) view = otv.View(graph) # %% # Other covariance models # ----------------------- # # There are other covariance models. The models which are used more often are the following: # # * :class:`~openturns.SquaredExponential`. The generated processes can be derivated in mean square at all orders. # * :class:`~openturns.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. # * :class:`~openturns.ExponentialModel`. The associated process is continuous, but not differentiable. # # %% # The Matérn and exponential models # --------------------------------- # %% 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) # %% nbTrajectories = 10 graph1 = plotCovarianceModel(myModel1, myTimeGrid, nbTrajectories) graph2 = plotCovarianceModel(myModel2, myTimeGrid, nbTrajectories) graph3 = plotCovarianceModel(myModel3, myTimeGrid, nbTrajectories) # %% fig = plt.figure(figsize=(20, 6)) ax1 = fig.add_subplot(1, 3, 1) _ = otv.View(graph1, figure=fig, axes=[ax1]) _ = ax1.set_title("Matern 5/2") ax2 = fig.add_subplot(1, 3, 2) _ = otv.View(graph2, figure=fig, axes=[ax2]) _ = ax2.set_title("Matern 3/2") ax3 = fig.add_subplot(1, 3, 3) _ = otv.View(graph3, figure=fig, axes=[ax3]) _ = ax3.set_title("Matern 1/2") # %% # We see than, when :math:`\nu` increases, then the trajectories are smoother and smoother. # %% myExpModel = ot.ExponentialModel(scale, amplitude) # %% graph = plotCovarianceModel(myExpModel, myTimeGrid, nbTrajectories) graph.setTitle("Exponential") view = otv.View(graph) # %% # We see that the exponential model produces very irregular trajectories. otv.View.ShowAll()