.. DO NOT EDIT. .. THIS FILE WAS AUTOMATICALLY GENERATED BY SPHINX-GALLERY. .. TO MAKE CHANGES, EDIT THE SOURCE PYTHON FILE: .. "auto_meta_modeling/kriging_metamodel/plot_kriging_sequential.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_meta_modeling_kriging_metamodel_plot_kriging_sequential.py: Sequentially adding new points to a kriging =========================================== .. GENERATED FROM PYTHON SOURCE LINES 6-8 In this example, we show how to sequentially add new points to a kriging in order to improve the predictivity of the metamodel. In order to create simple graphics, we consider a 1D function. .. GENERATED FROM PYTHON SOURCE LINES 10-12 Create the function and the design of experiments ------------------------------------------------- .. GENERATED FROM PYTHON SOURCE LINES 14-21 .. code-block:: Python import openturns as ot from openturns.viewer import View import numpy as np import openturns.viewer as viewer ot.Log.Show(ot.Log.NONE) .. GENERATED FROM PYTHON SOURCE LINES 22-25 .. code-block:: Python sampleSize = 4 dimension = 1 .. GENERATED FROM PYTHON SOURCE LINES 26-27 Define the function. .. GENERATED FROM PYTHON SOURCE LINES 29-31 .. code-block:: Python g = ot.SymbolicFunction(["x"], ["0.5*x^2 + sin(2.5*x)"]) .. GENERATED FROM PYTHON SOURCE LINES 32-33 Create the design of experiments. .. GENERATED FROM PYTHON SOURCE LINES 35-41 .. code-block:: Python xMin = -0.9 xMax = 1.9 X_distr = ot.Uniform(xMin, xMax) X = ot.LHSExperiment(X_distr, sampleSize, False, False).generate() Y = g(X) .. GENERATED FROM PYTHON SOURCE LINES 42-49 .. code-block:: Python graph = g.draw(xMin, xMax) data = ot.Cloud(X, Y) data.setColor("red") graph.add(data) view = viewer.View(graph) .. image-sg:: /auto_meta_modeling/kriging_metamodel/images/sphx_glr_plot_kriging_sequential_001.png :alt: y0 as a function of x :srcset: /auto_meta_modeling/kriging_metamodel/images/sphx_glr_plot_kriging_sequential_001.png :class: sphx-glr-single-img .. GENERATED FROM PYTHON SOURCE LINES 50-52 Create the algorithms --------------------- .. GENERATED FROM PYTHON SOURCE LINES 55-68 .. code-block:: Python def createMyBasicKriging(X, Y): """ Create a kriging from a pair of X and Y samples. We use a 3/2 Matérn covariance model and a constant trend. """ basis = ot.ConstantBasisFactory(dimension).build() covarianceModel = ot.MaternModel([1.0], 1.5) algo = ot.KrigingAlgorithm(X, Y, covarianceModel, basis) algo.run() krigResult = algo.getResult() return krigResult .. GENERATED FROM PYTHON SOURCE LINES 69-78 .. code-block:: Python def linearSample(xmin, xmax, npoints): """Returns a sample created from a regular grid from xmin to xmax with npoints points.""" step = (xmax - xmin) / (npoints - 1) rg = ot.RegularGrid(xmin, step, npoints) vertices = rg.getVertices() return vertices .. GENERATED FROM PYTHON SOURCE LINES 79-96 .. code-block:: Python def plot_kriging_bounds(vLow, vUp, n_test): """ From two lists containing the lower and upper bounds of the region, create a PolygonArray. """ palette = ot.Drawable.BuildDefaultPalette(2) myPaletteColor = palette[1] polyData = [[vLow[i], vLow[i + 1], vUp[i + 1], vUp[i]] for i in range(n_test - 1)] polygonList = [ ot.Polygon(polyData[i], myPaletteColor, myPaletteColor) for i in range(n_test - 1) ] boundsPoly = ot.PolygonArray(polygonList) boundsPoly.setLegend("95% bounds") return boundsPoly .. GENERATED FROM PYTHON SOURCE LINES 97-98 The following `sqrt` function will be used later to compute the standard deviation from the variance. .. GENERATED FROM PYTHON SOURCE LINES 100-103 .. code-block:: Python sqrt = ot.SymbolicFunction(["x"], ["sqrt(x)"]) .. GENERATED FROM PYTHON SOURCE LINES 104-168 .. code-block:: Python def plotMyBasicKriging(krigResult, xMin, xMax, X, Y, level=0.95): """ Given a kriging result, plot the data, the kriging metamodel and a confidence interval. """ samplesize = X.getSize() meta = krigResult.getMetaModel() graphKriging = meta.draw(xMin, xMax) graphKriging.setLegends(["Kriging"]) # Create a grid of points and evaluate the function and the kriging nbpoints = 50 xGrid = linearSample(xMin, xMax, nbpoints) yFunction = g(xGrid) yKrig = meta(xGrid) # Compute the conditional covariance epsilon = ot.Sample(nbpoints, [1.0e-8]) conditionalVariance = krigResult.getConditionalMarginalVariance(xGrid) + epsilon conditionalSigma = sqrt(conditionalVariance) # Compute the quantile of the Normal distribution alpha = 1 - (1 - level) / 2 quantileAlpha = ot.DistFunc.qNormal(alpha) # Graphics of the bounds epsilon = 1.0e-8 dataLower = [ yKrig[i, 0] - quantileAlpha * conditionalSigma[i, 0] for i in range(nbpoints) ] dataUpper = [ yKrig[i, 0] + quantileAlpha * conditionalSigma[i, 0] for i in range(nbpoints) ] # Coordinates of the vertices of the Polygons vLow = [[xGrid[i, 0], dataLower[i]] for i in range(nbpoints)] vUp = [[xGrid[i, 0], dataUpper[i]] for i in range(nbpoints)] # Compute the Polygon graphics boundsPoly = plot_kriging_bounds(vLow, vUp, nbpoints) boundsPoly.setLegend("95% bounds") # Validate the kriging metamodel mmv = ot.MetaModelValidation(xGrid, yFunction, meta) Q2 = mmv.computePredictivityFactor()[0] # Plot the function graphFonction = ot.Curve(xGrid, yFunction) graphFonction.setLineStyle("dashed") graphFonction.setColor("magenta") graphFonction.setLineWidth(2) graphFonction.setLegend("Function") # Draw the X and Y observed cloudDOE = ot.Cloud(X, Y) cloudDOE.setPointStyle("circle") cloudDOE.setColor("red") cloudDOE.setLegend("Data") # Assemble the graphics graph = ot.Graph() graph.add(boundsPoly) graph.add(graphFonction) graph.add(cloudDOE) graph.add(graphKriging) graph.setLegendPosition("lower right") graph.setAxes(True) graph.setGrid(True) graph.setTitle("Size = %d, Q2=%.2f%%" % (samplesize, 100 * Q2)) graph.setXTitle("X") graph.setYTitle("Y") return graph .. GENERATED FROM PYTHON SOURCE LINES 169-170 We start by creating the initial kriging metamodel on the 4 points in the design of experiments. .. GENERATED FROM PYTHON SOURCE LINES 172-177 .. code-block:: Python krigResult = createMyBasicKriging(X, Y) graph = plotMyBasicKriging(krigResult, xMin, xMax, X, Y) view = viewer.View(graph) .. image-sg:: /auto_meta_modeling/kriging_metamodel/images/sphx_glr_plot_kriging_sequential_002.png :alt: Size = 4, Q2=96.34% :srcset: /auto_meta_modeling/kriging_metamodel/images/sphx_glr_plot_kriging_sequential_002.png :class: sphx-glr-single-img .. GENERATED FROM PYTHON SOURCE LINES 178-180 Sequentially add new points --------------------------- .. GENERATED FROM PYTHON SOURCE LINES 182-183 The following function is the building block of the algorithm. It returns a new point which maximizes the conditional variance. .. GENERATED FROM PYTHON SOURCE LINES 186-200 .. code-block:: Python def getNewPoint(xMin, xMax, krigResult): """ Returns a new point to be added to the design of experiments. This point maximizes the conditional variance of the kriging. """ nbpoints = 50 xGrid = linearSample(xMin, xMax, nbpoints) conditionalVariance = krigResult.getConditionalMarginalVariance(xGrid) iMaxVar = int(np.argmax(conditionalVariance)) xNew = xGrid[iMaxVar, 0] xNew = ot.Point([xNew]) return xNew .. GENERATED FROM PYTHON SOURCE LINES 201-202 We first call `getNewPoint` to get a point to add to the design of experiments. .. GENERATED FROM PYTHON SOURCE LINES 205-208 .. code-block:: Python xNew = getNewPoint(xMin, xMax, krigResult) xNew .. raw:: html
class=Point name=Unnamed dimension=1 values=[1.9]


.. GENERATED FROM PYTHON SOURCE LINES 209-210 Then we evaluate the function on the new point and add it to the training design of experiments. .. GENERATED FROM PYTHON SOURCE LINES 212-216 .. code-block:: Python yNew = g(xNew) X.add(xNew) Y.add(yNew) .. GENERATED FROM PYTHON SOURCE LINES 217-218 We now plot the updated kriging. .. GENERATED FROM PYTHON SOURCE LINES 220-221 sphinx_gallery_thumbnail_number = 3 .. GENERATED FROM PYTHON SOURCE LINES 221-226 .. code-block:: Python krigResult = createMyBasicKriging(X, Y) graph = plotMyBasicKriging(krigResult, xMin, xMax, X, Y) graph.setTitle("Kriging #0") view = viewer.View(graph) .. image-sg:: /auto_meta_modeling/kriging_metamodel/images/sphx_glr_plot_kriging_sequential_003.png :alt: Kriging #0 :srcset: /auto_meta_modeling/kriging_metamodel/images/sphx_glr_plot_kriging_sequential_003.png :class: sphx-glr-single-img .. GENERATED FROM PYTHON SOURCE LINES 227-228 The algorithm added a point to the right bound of the domain. .. GENERATED FROM PYTHON SOURCE LINES 230-240 .. code-block:: Python for krigingStep in range(5): xNew = getNewPoint(xMin, xMax, krigResult) yNew = g(xNew) X.add(xNew) Y.add(yNew) krigResult = createMyBasicKriging(X, Y) graph = plotMyBasicKriging(krigResult, xMin, xMax, X, Y) graph.setTitle("Kriging #%d " % (krigingStep + 1) + graph.getTitle()) View(graph) .. rst-class:: sphx-glr-horizontal * .. image-sg:: /auto_meta_modeling/kriging_metamodel/images/sphx_glr_plot_kriging_sequential_004.png :alt: Kriging #1 Size = 6, Q2=98.80% :srcset: /auto_meta_modeling/kriging_metamodel/images/sphx_glr_plot_kriging_sequential_004.png :class: sphx-glr-multi-img * .. image-sg:: /auto_meta_modeling/kriging_metamodel/images/sphx_glr_plot_kriging_sequential_005.png :alt: Kriging #2 Size = 7, Q2=99.74% :srcset: /auto_meta_modeling/kriging_metamodel/images/sphx_glr_plot_kriging_sequential_005.png :class: sphx-glr-multi-img * .. image-sg:: /auto_meta_modeling/kriging_metamodel/images/sphx_glr_plot_kriging_sequential_006.png :alt: Kriging #3 Size = 8, Q2=99.87% :srcset: /auto_meta_modeling/kriging_metamodel/images/sphx_glr_plot_kriging_sequential_006.png :class: sphx-glr-multi-img * .. image-sg:: /auto_meta_modeling/kriging_metamodel/images/sphx_glr_plot_kriging_sequential_007.png :alt: Kriging #4 Size = 9, Q2=99.89% :srcset: /auto_meta_modeling/kriging_metamodel/images/sphx_glr_plot_kriging_sequential_007.png :class: sphx-glr-multi-img * .. image-sg:: /auto_meta_modeling/kriging_metamodel/images/sphx_glr_plot_kriging_sequential_008.png :alt: Kriging #5 Size = 10, Q2=99.90% :srcset: /auto_meta_modeling/kriging_metamodel/images/sphx_glr_plot_kriging_sequential_008.png :class: sphx-glr-multi-img .. GENERATED FROM PYTHON SOURCE LINES 241-245 We observe that the second added point is the left bound of the domain. The remaining points were added strictly inside the domain where the accuracy was drastically improved. With only 10 points, the metamodel accuracy is already very good with a Q2 which is equal to 99.9%. .. GENERATED FROM PYTHON SOURCE LINES 247-257 Conclusion ---------- The current example presents the naive implementation on the creation of a sequential design of experiments based on kriging. More practical algorithms are presented in the following references. * Mona Abtini. Plans prédictifs à taille fixe et séquentiels pour le krigeage (2008). Thèse de doctorat de l'Université de Lyon. * Céline Scheidt. Analyse statistique d’expériences simulées : Modélisation adaptative de réponses non régulières par krigeage et plans d’expériences (2007). Thèse présentée pour obtenir le grade de Docteur de l’Université Louis Pasteur. * David Ginsbourger. Sequential Design of Computer Experiments. Wiley StatsRef: Statistics Reference Online, Wiley (2018) .. GENERATED FROM PYTHON SOURCE LINES 257-259 .. code-block:: Python View.ShowAll() .. _sphx_glr_download_auto_meta_modeling_kriging_metamodel_plot_kriging_sequential.py: .. only:: html .. container:: sphx-glr-footer sphx-glr-footer-example .. container:: sphx-glr-download sphx-glr-download-jupyter :download:`Download Jupyter notebook: plot_kriging_sequential.ipynb ` .. container:: sphx-glr-download sphx-glr-download-python :download:`Download Python source code: plot_kriging_sequential.py `