.. 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_gpr_1d.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_gpr_1d.py: Gaussian Process Regression : quick-start ========================================= .. GENERATED FROM PYTHON SOURCE LINES 7-15 Abstract -------- In this example, we create a Gaussian process Regression for a function which has scalar real inputs and outputs. We show how to create the learning and the validation samples. We show how to create the metamodel by choosing a trend and a covariance model. Finally, we compute the predicted confidence interval using the conditional variance. .. GENERATED FROM PYTHON SOURCE LINES 17-48 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 Gaussian Process Regression metamodel with: * a constant trend, * a Matern covariance model. .. GENERATED FROM PYTHON SOURCE LINES 50-56 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. .. GENERATED FROM PYTHON SOURCE LINES 58-65 .. code-block:: Python import openturns as ot from openturns import viewer from matplotlib import pylab as plt import openturns.experimental as otexp ot.Log.Show(ot.Log.NONE) .. GENERATED FROM PYTHON SOURCE LINES 66-68 .. code-block:: Python g = ot.SymbolicFunction(["x"], ["sin(x)"]) .. GENERATED FROM PYTHON SOURCE LINES 69-74 .. code-block:: Python x_train = ot.Sample([[x] for x in [1.0, 3.0, 4.0, 6.0, 7.9, 11.0, 11.5]]) y_train = g(x_train) n_train = x_train.getSize() n_train .. rst-class:: sphx-glr-script-out .. code-block:: none 7 .. GENERATED FROM PYTHON SOURCE LINES 75-77 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. .. GENERATED FROM PYTHON SOURCE LINES 79-88 .. code-block:: Python xmin = 0.0 xmax = 12.0 n_test = 100 step = (xmax - xmin) / (n_test - 1) myRegularGrid = ot.RegularGrid(xmin, step, n_test) x_test = myRegularGrid.getVertices() y_test = g(x_test) .. GENERATED FROM PYTHON SOURCE LINES 89-90 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. .. GENERATED FROM PYTHON SOURCE LINES 92-108 .. code-block:: Python def plot_1d_data(x_data, y_data, type="Curve", legend=None, color=None, linestyle=None): """Plot the data (x_data,y_data) as a Cloud/Curve""" if type == "Curve": graphF = ot.Curve(x_data, y_data) else: graphF = ot.Cloud(x_data, y_data) if legend is not None: graphF.setLegend(legend) if color is not None: graphF.setColor(color) if linestyle is not None: graphF.setLineStyle(linestyle) return graphF .. GENERATED FROM PYTHON SOURCE LINES 109-118 .. code-block:: Python graph = ot.Graph("test and train", "", "", True, "") graph.add(plot_1d_data(x_test, y_test, legend="Exact", color="black", linestyle="dashed")) graph.add(plot_1d_data(x_train, y_train, type="Cloud", legend="Data", color="red")) graph.setAxes(True) graph.setXTitle("X") graph.setYTitle("Y") graph.setLegendPosition("upper right") view = viewer.View(graph) .. image-sg:: /auto_meta_modeling/kriging_metamodel/images/sphx_glr_plot_gpr_1d_001.png :alt: test and train :srcset: /auto_meta_modeling/kriging_metamodel/images/sphx_glr_plot_gpr_1d_001.png :class: sphx-glr-single-img .. GENERATED FROM PYTHON SOURCE LINES 119-121 We use the :class:`~openturns.ConstantBasisFactory` class to define the trend and the :class:`~openturns.MaternModel` class to define the covariance model. This Matérn model is based on the regularity parameter :math:`\nu=3/2`. .. GENERATED FROM PYTHON SOURCE LINES 123-134 .. code-block:: Python dimension = 1 basis = ot.ConstantBasisFactory(dimension).build() covarianceModel = ot.MaternModel([1.0] * dimension, 1.5) fitter_algo = otexp.GaussianProcessFitter(x_train, y_train, covarianceModel, basis) fitter_algo.run() fitter_result = fitter_algo.getResult() gpr_algo = otexp.GaussianProcessRegression(fitter_result) gpr_algo.run() gpr_result = gpr_algo.getResult() print(gpr_result) .. rst-class:: sphx-glr-script-out .. code-block:: none GaussianProcessRegressionResult(covariance models=MaternModel(scale=[1.27453], amplitude=[0.822263], nu=1.5), covariance coefficients=0 : [ 1.13904 ] 1 : [ 1.01762 ] 2 : [ -1.76279 ] 3 : [ -0.559148 ] 4 : [ 1.78757 ] 5 : [ -1.61946 ] 6 : [ -0.00283147 ], basis=Basis( [class=LinearEvaluation name=Unnamed center=[0] constant=[1] linear=[[ 0 ]]] ), trend coefficients=[0.00736753]) .. GENERATED FROM PYTHON SOURCE LINES 135-137 We observe that the `scale` and `amplitude` hyper-parameters have been optimized by the :meth:`~openturns.experimental.GaussianProcessFitter.run` method. Then we get the metamodel with `getMetaModel` and evaluate the outputs of the metamodel on the test design of experiments. .. GENERATED FROM PYTHON SOURCE LINES 139-143 .. code-block:: Python gprMetamodel = gpr_result.getMetaModel() y_test_MM = gprMetamodel(x_test) .. GENERATED FROM PYTHON SOURCE LINES 144-145 Now we plot Gaussian process Regression output, in addition to the previous plots .. GENERATED FROM PYTHON SOURCE LINES 148-158 .. code-block:: Python graph = ot.Graph("", "", "", True, "") graph.add(plot_1d_data(x_test, y_test, legend="Exact", color="black", linestyle="dashed")) graph.add(plot_1d_data(x_train, y_train, type="Cloud", legend="Data", color="red")) graph.add(plot_1d_data(x_test, y_test_MM, legend="GPR", color="blue")) graph.setAxes(True) graph.setXTitle("X") graph.setYTitle("Y") graph.setLegendPosition("upper right") view = viewer.View(graph) .. image-sg:: /auto_meta_modeling/kriging_metamodel/images/sphx_glr_plot_gpr_1d_002.png :alt: plot gpr 1d :srcset: /auto_meta_modeling/kriging_metamodel/images/sphx_glr_plot_gpr_1d_002.png :class: sphx-glr-single-img .. GENERATED FROM PYTHON SOURCE LINES 159-165 We see that the Gaussian process regression is interpolating. This is what is meant by *conditioning* a Gaussian process. We see that, when the sine function has a strong curvature between two points which are separated by a large distance (e.g. between :math:`x=4` and :math:`x=6`), then the gaussian regression is not close to the function :math:`g`. However, when the training points are close (e.g. between :math:`x=11` and :math:`x=11.5`) or when the function is nearly linear (e.g. between :math:`x=8` and :math:`x=11`), then the gaussian process regression is quite accurate. .. GENERATED FROM PYTHON SOURCE LINES 167-169 Compute confidence bounds ------------------------- .. GENERATED FROM PYTHON SOURCE LINES 171-180 In order to assess the quality of the metamodel, we can estimate the variance and compute a 95% confidence interval associated with the conditioned Gaussian process. We begin by defining the `alpha` variable containing the complementary of the confidence level than we want to compute. Then we compute the quantile of the Gaussian distribution corresponding to `1-alpha/2`. Therefore, the confidence interval is: .. math:: P\in\left(X\in\left[q_{\alpha/2},q_{1-\alpha/2}\right]\right)=1-\alpha. .. GENERATED FROM PYTHON SOURCE LINES 182-194 .. code-block:: Python alpha = 0.05 def computeQuantileAlpha(alpha): bilateralCI = ot.Normal().computeBilateralConfidenceInterval(1 - alpha) return bilateralCI.getUpperBound()[0] quantileAlpha = computeQuantileAlpha(alpha) print("alpha=%f" % (alpha)) print("Quantile alpha=%f" % (quantileAlpha)) .. rst-class:: sphx-glr-script-out .. code-block:: none alpha=0.050000 Quantile alpha=1.959964 .. GENERATED FROM PYTHON SOURCE LINES 195-204 In order to compute the regression error, we can consider the conditional variance. The :meth:`~openturns.experimental.GaussianProcessConditionalCovariance.getConditionalMarginalVariance` method returns the covariance matrix `covGrid` evaluated at each points in the given sample. Then we can use the diagonal coefficients in order to get the marginal conditional Kriging variance. Since this is a variance, we use the square root in order to compute the standard deviation. However, some coefficients in the diagonal are very close to zero and nonpositive, which leads to an exception of the sqrt function. This is why we add an epsilon on the diagonal (nugget factor), which prevents this issue. .. GENERATED FROM PYTHON SOURCE LINES 206-212 .. code-block:: Python sqrt = ot.SymbolicFunction(["x"], ["sqrt(x)"]) epsilon = ot.Sample(n_test, [1.0e-8]) gccc = otexp.GaussianProcessConditionalCovariance(gpr_result) conditionalVariance = gccc.getConditionalMarginalVariance(x_test) + epsilon conditionalSigma = sqrt(conditionalVariance) .. GENERATED FROM PYTHON SOURCE LINES 213-214 The following figure presents the conditional standard deviation depending on :math:`x`. .. GENERATED FROM PYTHON SOURCE LINES 216-224 .. code-block:: Python graph = ot.Graph( "Conditional standard deviation", "x", "Conditional standard deviation", True, "" ) curve = ot.Curve(x_test, conditionalSigma) graph.add(curve) view = viewer.View(graph) .. image-sg:: /auto_meta_modeling/kriging_metamodel/images/sphx_glr_plot_gpr_1d_003.png :alt: Conditional standard deviation :srcset: /auto_meta_modeling/kriging_metamodel/images/sphx_glr_plot_gpr_1d_003.png :class: sphx-glr-single-img .. GENERATED FROM PYTHON SOURCE LINES 225-227 We now compute the bounds of the confidence interval. For this purpose we define a small function `computeBoundsConfidenceInterval` : .. GENERATED FROM PYTHON SOURCE LINES 230-244 .. code-block:: Python def computeBoundsConfidenceInterval(quantileAlpha): dataLower = [ [y_test_MM[i, 0] - quantileAlpha * conditionalSigma[i, 0]] for i in range(n_test) ] dataUpper = [ [y_test_MM[i, 0] + quantileAlpha * conditionalSigma[i, 0]] for i in range(n_test) ] dataLower = ot.Sample(dataLower) dataUpper = ot.Sample(dataUpper) return dataLower, dataUpper .. GENERATED FROM PYTHON SOURCE LINES 245-246 We define two small lists to draw three different confidence intervals (defined by the alpha value) : .. GENERATED FROM PYTHON SOURCE LINES 246-250 .. code-block:: Python alphas = [0.05, 0.1, 0.2] # three different green colors defined by HSV values mycolors = [[120, 1.0, 1.0], [120, 1.0, 0.75], [120, 1.0, 0.5]] .. GENERATED FROM PYTHON SOURCE LINES 251-252 We are ready to display all the previous information and the three confidence intervals we want. .. GENERATED FROM PYTHON SOURCE LINES 254-255 sphinx_gallery_thumbnail_number = 4 .. GENERATED FROM PYTHON SOURCE LINES 255-278 .. code-block:: Python graph = ot.Graph("", "", "", True, "") # Now we loop over the different values : for idx, v in enumerate(alphas): quantileAlpha = computeQuantileAlpha(v) vLow, vUp = computeBoundsConfidenceInterval(quantileAlpha) boundsPoly = ot.Polygon.FillBetween(x_test, vLow, vUp) boundsPoly.setColor( ot.Drawable.ConvertFromHSV(mycolors[idx][0], mycolors[idx][1], mycolors[idx][2]) ) boundsPoly.setLegend(" %d%% bounds" % ((1.0 - v) * 100)) graph.add(boundsPoly) graph.add(plot_1d_data(x_test, y_test, legend="Exact", color="black", linestyle="dashed")) graph.add(plot_1d_data(x_train, y_train, type="Cloud", legend="Data", color="red")) graph.add(plot_1d_data(x_test, y_test_MM, legend="GPR", color="blue")) graph.setAxes(True) graph.setXTitle("X") graph.setYTitle("Y") graph.setLegendPosition("upper right") view = viewer.View(graph) .. image-sg:: /auto_meta_modeling/kriging_metamodel/images/sphx_glr_plot_gpr_1d_004.png :alt: plot gpr 1d :srcset: /auto_meta_modeling/kriging_metamodel/images/sphx_glr_plot_gpr_1d_004.png :class: sphx-glr-single-img .. GENERATED FROM PYTHON SOURCE LINES 279-284 We see that the confidence intervals are small in the regions where two consecutive training points are close to each other (e.g. between :math:`x=11` and :math:`x=11.5`) and large when the two points are not (e.g. between :math:`x=8.` and :math:`x=11`) or when the curvature of the function is large (between :math:`x=4` and :math:`x=6`). .. GENERATED FROM PYTHON SOURCE LINES 284-286 .. code-block:: Python plt.show() .. _sphx_glr_download_auto_meta_modeling_kriging_metamodel_plot_gpr_1d.py: .. only:: html .. container:: sphx-glr-footer sphx-glr-footer-example .. container:: sphx-glr-download sphx-glr-download-jupyter :download:`Download Jupyter notebook: plot_gpr_1d.ipynb ` .. container:: sphx-glr-download sphx-glr-download-python :download:`Download Python source code: plot_gpr_1d.py ` .. container:: sphx-glr-download sphx-glr-download-zip :download:`Download zipped: plot_gpr_1d.zip `