.. only:: html .. note:: :class: sphx-glr-download-link-note Click :ref:`here ` to download the full example code .. rst-class:: sphx-glr-example-title .. _sphx_glr_auto_meta_modeling_kriging_metamodel_plot_kriging_simulate.py: Simulate new trajectories from a kriging metamodel ================================================== The main goal of this example is to show how to simulate new trajectories from a kriging metamodel. Introduction ------------ We consider the sine function: .. math:: y = \sin(x) for any :math:`x\in[0,12]`. We want to create a metamodel of this function. This is why we create a sample of :math:`n` observations of the function: .. math:: y_i=\sin(x_i) for :math:`i=1,...,7`, where :math:`x_i` is the i-th input and :math:`y_i` is the corresponding output. We consider the seven following inputs : ============ === === === === ===== ==== ====== :math:`i` 1 2 3 4 5 6 7 :math:`x_i` 1 3 4 6 7.9 11 11.5 ============ === === === === ===== ==== ====== We are going to consider a kriging metamodel with a * constant trend, * a Matern covariance model. Creation of the metamodel ------------------------- We begin by defining the function `g` as a symbolic function. Then we define the `x_train` variable which contains the inputs of the design of experiments of the training step. Then we compute the `y_train` corresponding outputs. The variable `n_train` is the size of the training design of experiments. .. code-block:: default import numpy as np import openturns as ot import openturns.viewer as viewer from matplotlib import pylab as plt ot.Log.Show(ot.Log.NONE) .. code-block:: default g = ot.SymbolicFunction(['x'], ['sin(x)']) .. code-block:: default x_train = ot.Sample([[x] for x in [1.,3.,4.,6.,7.9,11., 11.5]]) y_train = g(x_train) n_train = x_train.getSize() n_train .. rst-class:: sphx-glr-script-out Out: .. code-block:: none 7 In order to compare the function and its metamodel, we use a test (i.e. validation) design of experiments made of a regular grid of 100 points from 0 to 12. Then we convert this grid into a `Sample` and we compute the outputs of the function on this sample. .. code-block:: default xmin = 0. xmax = 12. n_test = 100 step = (xmax-xmin)/(n_test-1) myRegularGrid = ot.RegularGrid(xmin, step, n_test) x_test_coord = myRegularGrid.getValues() x_test = ot.Sample([[x] for x in x_test_coord]) y_test = g(x_test) In order to observe the function and the location of the points in the input design of experiments, we define the following functions which plots the data. .. code-block:: default def plot_data_train(x_train,y_train): '''Plot the data (x_train,y_train) as a Cloud, in red''' graph_train = ot.Cloud(x_train,y_train) graph_train.setColor("red") graph_train.setLegend("Data") return graph_train .. code-block:: default def plot_data_test(x_test,y_test): '''Plot the data (x_test,y_test) as a Curve, in dashed black''' graphF = ot.Curve(x_test,y_test) graphF.setLegend("Exact") graphF.setColor("black") graphF.setLineStyle("dashed") return graphF .. code-block:: default graph = ot.Graph() graph.add(plot_data_test(x_test,y_test)) graph.add(plot_data_train(x_train,y_train)) graph.setAxes(True) graph.setXTitle("X") graph.setYTitle("Y") graph.setLegendPosition("topright") view = viewer.View(graph) .. image:: /auto_meta_modeling/kriging_metamodel/images/sphx_glr_plot_kriging_simulate_001.png :alt: plot kriging simulate :class: sphx-glr-single-img We use the `ConstantBasisFactory` class to define the trend and the `MaternModel` class to define the covariance model. This Matérn model is based on the regularity parameter :math:`\nu=3/2`. .. code-block:: default dimension = 1 basis = ot.ConstantBasisFactory(dimension).build() covarianceModel = ot.MaternModel([1.]*dimension, 1.5) algo = ot.KrigingAlgorithm(x_train, y_train, covarianceModel, basis) algo.run() krigingResult = algo.getResult() krigingResult .. raw:: html

KrigingResult(covariance models=MaternModel(scale=[1.26541], amplitude=[0.818994], nu=1.5), covariance coefficients=0 : [ 1.14764 ]
1 : [ 1.006 ]
2 : [ -1.75286 ]
3 : [ -0.55873 ]
4 : [ 1.78675 ]
5 : [ -1.59663 ]
6 : [ -0.032173 ], basis=[Basis( [class=LinearEvaluation name=Unnamed center=[0] constant=[1] linear=[[ 0 ]]] )], trend coefficients=[[0.00671564]])



We observe that the `scale` and `amplitude` hyper-parameters have been optimized by the `run` method. Then we get the metamodel with `getMetaModel` and evaluate the outputs of the metamodel on the test design of experiments. .. code-block:: default krigeageMM = krigingResult.getMetaModel() y_test_MM = krigeageMM(x_test) The following function plots the kriging data. .. code-block:: default def plot_data_kriging(x_test,y_test_MM): '''Plots (x_test,y_test_MM) from the metamodel as a Curve, in blue''' graphK = ot.Curve(x_test,y_test_MM) graphK.setColor("blue") graphK.setLegend("Kriging") return graphK .. code-block:: default graph = ot.Graph() graph.add(plot_data_test(x_test,y_test)) graph.add(plot_data_train(x_train,y_train)) graph.add(plot_data_kriging(x_test,y_test_MM)) graph.setAxes(True) graph.setXTitle("X") graph.setYTitle("Y") graph.setLegendPosition("topright") view = viewer.View(graph) .. image:: /auto_meta_modeling/kriging_metamodel/images/sphx_glr_plot_kriging_simulate_002.png :alt: plot kriging simulate :class: sphx-glr-single-img Simulate new trajectories ------------------------- In order to generate new trajectories of the conditioned gaussian process, we couild technically use the `KrigingRandomVector` class, because it provides the `getSample` method that we need. However, the `KrigingRandomVector` class was more specifically designed to create a `RandomVector` so that it can feed, for example, a function which has a field as input argument. This is why we use the `ConditionedGaussianProcess`, which provides a `Process`. .. code-block:: default n_test = 100 step = (xmax-xmin)/(n_test-1) myRegularGrid = ot.RegularGrid(xmin, step, n_test) vertices = myRegularGrid.getVertices() If we directly use the `vertices` values, we get: RuntimeError: InternalException : Error: the matrix is not definite positive. Indeed, some points in `vertices` are also in `x_train`. This is why the conditioned covariance matrix is singular at these points. This is why we define the following function which deletes points in `vertices` which are also found in `x_train`. .. code-block:: default def deleteCommonValues(x_train,x_test): ''' Delete from x_test the values which are in x_train so that values in x_test have no interect with x_train. ''' x_test_filtered = x_test # Initialize for x_train_value in x_train: print("Checking %s" % (x_train_value)) indices = np.argwhere(x_test==x_train_value) if len(indices) == 1: print(" Delete %s" % (x_train_value)) x_test_filtered = np.delete(x_test_filtered, indices[0, 0]) else: print(" OK") return x_test_filtered .. code-block:: default vertices_filtered = deleteCommonValues(np.array(x_train.asPoint()),np.array(vertices.asPoint())) .. rst-class:: sphx-glr-script-out Out: .. code-block:: none Checking 1.0 OK Checking 3.0 OK Checking 4.0 Delete 4.0 Checking 6.0 OK Checking 7.9 OK Checking 11.0 OK Checking 11.5 OK .. code-block:: default evaluationMesh = ot.Mesh(ot.Sample([[vf] for vf in vertices_filtered])) .. code-block:: default process = ot.ConditionedGaussianProcess(krigingResult, evaluationMesh) .. code-block:: default trajectories = process.getSample(10) type(trajectories) The `getSample` method returns a `ProcessSample`. By comparison, the `getSample` method of a `KrigingRandomVector` would return a `Sample`. .. code-block:: default graph = trajectories.drawMarginal() graph.add(plot_data_test(x_test,y_test)) graph.add(plot_data_train(x_train,y_train)) graph.setAxes(True) graph.setXTitle("X") graph.setYTitle("Y") graph.setLegendPosition("topright") graph.setTitle("10 simulated trajectories") view = viewer.View(graph) plt.show() .. image:: /auto_meta_modeling/kriging_metamodel/images/sphx_glr_plot_kriging_simulate_003.png :alt: 10 simulated trajectories :class: sphx-glr-single-img References ---------- * Metamodeling with Gaussian processes, Bertrand Iooss, EDF R&D, 2014, www.gdr-mascotnum.fr/media/sssamo14_iooss.pdf .. rst-class:: sphx-glr-timing **Total running time of the script:** ( 0 minutes 0.265 seconds) .. _sphx_glr_download_auto_meta_modeling_kriging_metamodel_plot_kriging_simulate.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_kriging_simulate.py ` .. container:: sphx-glr-download sphx-glr-download-jupyter :download:`Download Jupyter notebook: plot_kriging_simulate.ipynb ` .. only:: html .. rst-class:: sphx-glr-signature `Gallery generated by Sphinx-Gallery `_