.. DO NOT EDIT. .. THIS FILE WAS AUTOMATICALLY GENERATED BY SPHINX-GALLERY. .. TO MAKE CHANGES, EDIT THE SOURCE PYTHON FILE: .. "auto_calibration/least_squares_and_gaussian_calibration/plot_calibration_flooding.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_calibration_least_squares_and_gaussian_calibration_plot_calibration_flooding.py: Calibration of the flooding model ================================= In this example we are interested in the calibration of the :ref:`flooding model `. We calibrate the parameters of a flooding model where only the difference between the downstream and upstream riverbed levels can be calibrated. This example shows how to manage the lack of identifiability in a calibration problem. This example use least squares to calibrate the parametric model. Please read :ref:`code_calibration` for more details on code calibration and least squares. This study is relatively complicated: please read the :doc:`calibration of the Chaboche mechanical model ` first if this is not already done. The observations that we use in this study are simulated with the script :doc:`Generate flooding model observations `. .. GENERATED FROM PYTHON SOURCE LINES 20-48 Parameters to calibrate and observations ---------------------------------------- The variables of the model are: - :math:`Q` : Input. Observed. - :math:`K_s`, :math:`Z_v`, :math:`Z_m` : Input. Calibrated. - :math:`H`: Output. Observed. The vector of parameters to calibrate is: .. math:: \theta = (K_s,Z_v,Z_m). In the description of the :ref:`flooding model`, we see that only one parameter can be identified. Hence, calibrating this model requires some regularization. We return to this topic when analyzing the singular values of the Jacobian matrix. We consider a sample size equal to: .. math:: n = 10. The observations are the couples :math:`\{(Q_i,H_i)\}_{i=1,...,n}`, i.e. each observation is a couple made of the flowrate and the corresponding river height. .. GENERATED FROM PYTHON SOURCE LINES 48-58 .. code-block:: Python from openturns.usecases import flood_model from matplotlib import pylab as plt import openturns.viewer as otv import numpy as np import openturns as ot ot.ResourceMap.SetAsUnsignedInteger("Normal-SmallDimension", 1) ot.Log.Show(ot.Log.NONE) .. GENERATED FROM PYTHON SOURCE LINES 59-67 Define the observations ----------------------- In practice, we generally use a data set which has been obtained from measurements. This data set can be loaded using e.g. :meth:`~openturns.Sample.ImportFromCSVFile`. Here we import the data from the :class:`~openturns.usecases.flood_model.FloodModel` class. .. GENERATED FROM PYTHON SOURCE LINES 67-73 .. code-block:: Python fm = flood_model.FloodModel() Qobs = fm.data[:, 0] Hobs = fm.data[:, 1] nbobs = fm.data.getSize() fm.data .. raw:: html
Q ($m^3/s$)H (m)
01300.59
15301.33
29602.03
314002.72
418302.83
522603.5
627003.82
731304.36
835604.63
940104.96


.. GENERATED FROM PYTHON SOURCE LINES 74-88 Create the physical model ------------------------- We define the model :math:`g` which has 4 inputs and one output H. The nonlinear least squares algorithm does not take into account for bounds in the parameters. Therefore, we ensure that the output is computed whatever the inputs. The model fails into two situations: * if :math:`K_s<0`, * if :math:`Z_v-Z_m<0`. In these cases, we return an infinite number, so that the optimization algorithm does not get trapped. .. GENERATED FROM PYTHON SOURCE LINES 90-109 .. code-block:: Python def functionFlooding(X): L = 5.0e3 B = 300.0 Q, K_s, Z_v, Z_m = X alpha = (Z_m - Z_v) / L if alpha < 0.0 or K_s <= 0.0: H = np.inf else: H = (Q / (K_s * B * np.sqrt(alpha))) ** (3.0 / 5.0) return [H] g = ot.PythonFunction(4, 1, functionFlooding) g = ot.MemoizeFunction(g) g.setInputDescription(["Q ($m^3/s$)", "Ks ($m^{1/3}/s$)", "Zv (m)", "Zm (m)"]) g.setOutputDescription(["H (m)"]) .. GENERATED FROM PYTHON SOURCE LINES 110-116 Setting the calibration parameters ---------------------------------- Define the value of the reference values of the :math:`\theta` parameter. In the Bayesian framework, this is called the mean of the *prior* normal distribution. In the data assimilation framework, this is called the *background*. .. GENERATED FROM PYTHON SOURCE LINES 116-121 .. code-block:: Python KsInitial = 20.0 ZvInitial = 49.0 ZmInitial = 51.0 thetaPrior = [KsInitial, ZvInitial, ZmInitial] .. GENERATED FROM PYTHON SOURCE LINES 122-148 Create the parametric function ------------------------------ In the physical model, the inputs and parameters are ordered as presented in the next table. Notice that there are no parameters in the physical model. +-------+----------------+ | Index | Input variable | +=======+================+ | 0 | Q | +-------+----------------+ | 1 | Ks | +-------+----------------+ | 2 | Zv | +-------+----------------+ | 3 | Zm | +-------+----------------+ +-------+-----------+ | Index | Parameter | +=======+===========+ | ∅ | ∅ | +-------+-----------+ **Table 1.** Indices and names of the inputs and parameters of the physical model. .. GENERATED FROM PYTHON SOURCE LINES 148-151 .. code-block:: Python print("Physical Model Inputs:", g.getInputDescription()) print("Physical Model Parameters:", g.getParameterDescription()) .. rst-class:: sphx-glr-script-out .. code-block:: none Physical Model Inputs: [Q ($m^3/s$),Ks ($m^{1/3}/s$),Zv (m),Zm (m)] Physical Model Parameters: [] .. GENERATED FROM PYTHON SOURCE LINES 152-176 In order to perform calibration, we have to define a parametric model, with observed inputs and parameters to calibrate. In order to do this, we create a :class:`~openturns.ParametricFunction` where the parameters are `Ks`, `Zv` and `Zm` which have the indices 1, 2 and 3 in the physical model. +-------+----------------+ | Index | Input variable | +=======+================+ | 0 | Q | +-------+----------------+ +-------+-----------+ | Index | Parameter | +=======+===========+ | 0 | Ks | +-------+-----------+ | 1 | Zv | +-------+-----------+ | 3 | Zm | +-------+-----------+ **Table 2.** Indices and names of the inputs and parameters of the parametric model. .. GENERATED FROM PYTHON SOURCE LINES 179-182 The following statement creates the calibrated function from the model. The calibrated parameters :math:`K_s`, :math:`Z_v`, :math:`Z_m` are at indices 1, 2, 3 in the inputs arguments of the model. .. GENERATED FROM PYTHON SOURCE LINES 182-185 .. code-block:: Python calibratedIndices = [1, 2, 3] mycf = ot.ParametricFunction(g, calibratedIndices, thetaPrior) .. GENERATED FROM PYTHON SOURCE LINES 186-187 Plot the Y observations versus the X observations. .. GENERATED FROM PYTHON SOURCE LINES 187-204 .. code-block:: Python graph = ot.Graph("Observations", "Q ($m^3/s$)", "H (m)", True) # Plot the model before calibration curve = mycf.draw(100.0, 4000.0).getDrawable(0) curve.setLegend("Model, before calibration") curve.setLineStyle(ot.ResourceMap.GetAsString("CalibrationResult-ObservationLineStyle")) graph.add(curve) # Plot the noisy observations cloud = ot.Cloud(Qobs, Hobs) cloud.setLegend("Observations") cloud.setPointStyle( ot.ResourceMap.GetAsString("CalibrationResult-ObservationPointStyle") ) graph.add(cloud) # graph.setLegendPosition("upper left") view = otv.View(graph) .. image-sg:: /auto_calibration/least_squares_and_gaussian_calibration/images/sphx_glr_plot_calibration_flooding_001.png :alt: Observations :srcset: /auto_calibration/least_squares_and_gaussian_calibration/images/sphx_glr_plot_calibration_flooding_001.png :class: sphx-glr-single-img .. GENERATED FROM PYTHON SOURCE LINES 205-208 Wee see that the model does not fit to the data. The goal of calibration is to find which parameter best fit to the observations. .. GENERATED FROM PYTHON SOURCE LINES 210-215 Calibration with linear least squares ------------------------------------- The :class:`~openturns.LinearLeastSquaresCalibration` class performs the linear least squares calibration by linearizing the model in the neighbourhood of the reference point. .. GENERATED FROM PYTHON SOURCE LINES 215-217 .. code-block:: Python algo = ot.LinearLeastSquaresCalibration(mycf, Qobs, Hobs, thetaPrior, "SVD") .. GENERATED FROM PYTHON SOURCE LINES 218-220 The :meth:`~openturns.LinearLeastSquaresCalibration.run` method computes the solution of the problem. .. GENERATED FROM PYTHON SOURCE LINES 220-223 .. code-block:: Python algo.run() calibrationResult = algo.getResult() .. GENERATED FROM PYTHON SOURCE LINES 224-226 The :meth:`~openturns.CalibrationResult.getParameterMAP` method returns the maximum of the posterior distribution of :math:`\theta`. .. GENERATED FROM PYTHON SOURCE LINES 226-229 .. code-block:: Python thetaMAP = calibrationResult.getParameterMAP() print("theta After = ", thetaMAP) print("theta Before = ", thetaPrior) .. rst-class:: sphx-glr-script-out .. code-block:: none theta After = [2.94129e+09,-3.26646e+24,-3.26646e+24] theta Before = [20.0, 49.0, 51.0] .. GENERATED FROM PYTHON SOURCE LINES 230-231 Print the true values of the parameters. .. GENERATED FROM PYTHON SOURCE LINES 231-236 .. code-block:: Python print("True theta") print(" Ks = ", fm.trueKs) print(" Zv = ", fm.trueZv) print(" Zm = ", fm.trueZm) .. rst-class:: sphx-glr-script-out .. code-block:: none True theta Ks = 30.0 Zv = 50.0 Zm = 55.0 .. GENERATED FROM PYTHON SOURCE LINES 237-241 In this case, we see that there seems to be a great distance from the reference value of :math:`\theta` to the optimum: the values seem too large in magnitude. As we are going to see, there is an identification problem because the Jacobian matrix is rank-degenerate. .. GENERATED FROM PYTHON SOURCE LINES 243-245 Diagnostic of the identification issue -------------------------------------- .. GENERATED FROM PYTHON SOURCE LINES 247-251 In this section, we show how to diagnose the identification problem. The :meth:`~openturns.CalibrationResult.getParameterPosterior` method returns the posterior normal distribution of :math:`\theta`. .. GENERATED FROM PYTHON SOURCE LINES 253-256 .. code-block:: Python distributionPosterior = calibrationResult.getParameterPosterior() print(distributionPosterior) .. rst-class:: sphx-glr-script-out .. code-block:: none Normal(mu = [2.94129e+09,-3.26646e+24,-3.26646e+24], sigma = [3.03385e+28,2.38506e+34,2.38506e+34], R = [[ 1 4.9106e-26 -4.9106e-26 ] [ 4.9106e-26 1 1 ] [ -4.9106e-26 1 1 ]]) .. GENERATED FROM PYTHON SOURCE LINES 257-260 We see that there is a large covariance matrix diagonal. Let us compute a 95% confidence interval for the solution :math:`\theta^\star`. .. GENERATED FROM PYTHON SOURCE LINES 262-268 .. code-block:: Python print( distributionPosterior.computeBilateralConfidenceIntervalWithMarginalProbability( 0.95 )[0] ) .. rst-class:: sphx-glr-script-out .. code-block:: none [-6.78405e+28, 6.78405e+28] [-5.33329e+34, 5.33329e+34] [-5.33329e+34, 5.33329e+34] .. GENERATED FROM PYTHON SOURCE LINES 269-272 The confidence interval is *very* large. In order to clarify the situation, we compute the Jacobian matrix of the model at the candidate point. .. GENERATED FROM PYTHON SOURCE LINES 274-281 .. code-block:: Python mycf.setParameter(thetaPrior) thetaDim = len(thetaPrior) jacobianMatrix = ot.Matrix(nbobs, thetaDim) for i in range(nbobs): jacobianMatrix[i, :] = mycf.parameterGradient(Qobs[i]).transpose() print(jacobianMatrix[0:5, :]) .. rst-class:: sphx-glr-script-out .. code-block:: none 5x3 [[ -0.0314759 0.15738 -0.15738 ] [ -0.0731439 0.365719 -0.365719 ] [ -0.104466 0.52233 -0.52233 ] [ -0.131005 0.655027 -0.655027 ] [ -0.153845 0.769225 -0.769225 ]] .. GENERATED FROM PYTHON SOURCE LINES 282-284 The rank of the problem can be seen from the singular values of the Jacobian matrix. .. GENERATED FROM PYTHON SOURCE LINES 286-288 .. code-block:: Python print(jacobianMatrix.computeSingularValues()) .. rst-class:: sphx-glr-script-out .. code-block:: none [3.8102,5.32438e-11,3.00517e-26] .. GENERATED FROM PYTHON SOURCE LINES 289-306 We can see that there are two singular values which are relatively close to zero. This explains why the Jacobian matrix is close to being rank-degenerate. Moreover, this allows one to compute the actual dimensionality of the problem. The algorithm we use computes the singular values in descending order. Moreover, by definition, the singular values are nonnegative. We see that the first singular value is close to :math:`10` and the others are very close to :math:`0` in comparison. This implies that the (numerical) rank of the Jacobian matrix is 1, even if there are 3 parameters. Hence, only one parameter can be identified, be it :math:`K_s`, :math:`Z_v` or :math:`Z_m`. The choice of the particular parameter to identify is free. However, in hydraulic studies, the parameter :math:`K_s` is classically calibrated while :math:`Z_v` and :math:`Z_m` are left constant. .. GENERATED FROM PYTHON SOURCE LINES 308-310 Conclusion of the linear least squares calibration -------------------------------------------------- .. GENERATED FROM PYTHON SOURCE LINES 312-324 There are several methods to solve the problem. * Given that the problem is not identifiable, we can use some regularization method. Two methods are provided in the library: the Gaussian linear least squares `GaussianLinearCalibration` and the Gaussian non linear least squares `GaussianNonlinearCalibration`. * We can change the problem, replacing it with a problem which is identifiable. In the flooding model, we can view :math:`Z_v` and :math:`Z_m` as constants and calibrate :math:`K_s` only. In this example, we do not change the problem and see how the different methods perform. .. GENERATED FROM PYTHON SOURCE LINES 326-331 Calibration with non linear least squares ----------------------------------------- The :class:`~openturns.NonLinearLeastSquaresCalibration` class performs the non linear least squares calibration by minimizing the squared Euclidian norm between the predictions and the observations. .. GENERATED FROM PYTHON SOURCE LINES 333-335 .. code-block:: Python algo = ot.NonLinearLeastSquaresCalibration(mycf, Qobs, Hobs, thetaPrior) .. GENERATED FROM PYTHON SOURCE LINES 336-338 The :meth:`~openturns.NonLinearLeastSquaresCalibration.run` method computes the solution of the problem. .. GENERATED FROM PYTHON SOURCE LINES 338-341 .. code-block:: Python algo.run() calibrationResult = algo.getResult() .. GENERATED FROM PYTHON SOURCE LINES 342-344 Analysis of the results ----------------------- .. GENERATED FROM PYTHON SOURCE LINES 346-348 The :meth:`~openturns.CalibrationResult.getParameterMAP` method returns the maximum of the posterior distribution of :math:`\theta`. .. GENERATED FROM PYTHON SOURCE LINES 350-353 .. code-block:: Python thetaMAP = calibrationResult.getParameterMAP() print(thetaMAP) .. rst-class:: sphx-glr-script-out .. code-block:: none [27.566,47.0918,52.9082] .. GENERATED FROM PYTHON SOURCE LINES 354-361 We can compute a 95% confidence interval of the parameter :math:`\theta^\star`. This confidence interval is based on bootstrap, based on a sample size equal to 100 (as long as the value of the :class:`~openturns.ResourceMap` key "NonLinearLeastSquaresCalibration-BootstrapSize" is unchanged). This confidence interval reflects the sensitivity of the optimum to the variability in the observations. .. GENERATED FROM PYTHON SOURCE LINES 361-367 .. code-block:: Python print( ot.ResourceMap.GetAsUnsignedInteger( "NonLinearLeastSquaresCalibration-BootstrapSize" ) ) .. rst-class:: sphx-glr-script-out .. code-block:: none 100 .. GENERATED FROM PYTHON SOURCE LINES 368-371 .. code-block:: Python thetaPosterior = calibrationResult.getParameterPosterior() print(thetaPosterior.computeBilateralConfidenceIntervalWithMarginalProbability(0.95)[0]) .. rst-class:: sphx-glr-script-out .. code-block:: none [27.2377, 28.1071] [46.863, 47.2304] [52.7696, 53.137] .. GENERATED FROM PYTHON SOURCE LINES 372-374 In this case, the values of the parameters are quite accurately computed. .. GENERATED FROM PYTHON SOURCE LINES 376-378 Increase the default number of points in the plots. This produces smoother spiky distributions. .. GENERATED FROM PYTHON SOURCE LINES 378-380 .. code-block:: Python ot.ResourceMap.SetAsUnsignedInteger("Distribution-DefaultPointNumber", 300) .. GENERATED FROM PYTHON SOURCE LINES 381-385 .. code-block:: Python graph = calibrationResult.drawObservationsVsInputs() graph.setLegendPosition("upper left") view = otv.View(graph) .. image-sg:: /auto_calibration/least_squares_and_gaussian_calibration/images/sphx_glr_plot_calibration_flooding_002.png :alt: plot calibration flooding :srcset: /auto_calibration/least_squares_and_gaussian_calibration/images/sphx_glr_plot_calibration_flooding_002.png :class: sphx-glr-single-img .. GENERATED FROM PYTHON SOURCE LINES 386-388 We see that there is a good fit after calibration, since the predictions after calibration are close to the observations. .. GENERATED FROM PYTHON SOURCE LINES 390-393 .. code-block:: Python graph = calibrationResult.drawObservationsVsPredictions() view = otv.View(graph) .. image-sg:: /auto_calibration/least_squares_and_gaussian_calibration/images/sphx_glr_plot_calibration_flooding_003.png :alt: plot calibration flooding :srcset: /auto_calibration/least_squares_and_gaussian_calibration/images/sphx_glr_plot_calibration_flooding_003.png :class: sphx-glr-single-img .. GENERATED FROM PYTHON SOURCE LINES 394-396 We see that there is a much better fit after calibration, since the predictions are close to the diagonal of the graphics. .. GENERATED FROM PYTHON SOURCE LINES 398-401 .. code-block:: Python observationError = calibrationResult.getObservationsError() print(observationError) .. rst-class:: sphx-glr-script-out .. code-block:: none Normal(mu = 0.0185511, sigma = 0.110646) .. GENERATED FROM PYTHON SOURCE LINES 402-404 We can see that the observation error has a sample mean close to zero and a sample standard deviation approximately equal to 0.11. .. GENERATED FROM PYTHON SOURCE LINES 404-410 .. code-block:: Python # sphinx_gallery_thumbnail_number = 5 graph = calibrationResult.drawResiduals() graph.setLegendPosition("upper left") view = otv.View(graph) .. image-sg:: /auto_calibration/least_squares_and_gaussian_calibration/images/sphx_glr_plot_calibration_flooding_004.png :alt: Residual analysis :srcset: /auto_calibration/least_squares_and_gaussian_calibration/images/sphx_glr_plot_calibration_flooding_004.png :class: sphx-glr-single-img .. GENERATED FROM PYTHON SOURCE LINES 411-418 The analysis of the residuals shows that the distribution is centered on zero and symmetric. This indicates that the calibration performed well. Moreover, the distribution of the residuals is close to being Gaussian. This is an important hypothesis of the least squares method so that checking that this hypothesis occurs in the study is an important verification. .. GENERATED FROM PYTHON SOURCE LINES 420-429 .. code-block:: Python graph = calibrationResult.drawParameterDistributions() view = otv.View( graph, figure_kw={"figsize": (8.0, 3.0)}, legend_kw={"bbox_to_anchor": (1.0, 1.0), "loc": "upper left"}, ) plt.subplots_adjust(right=0.8, bottom=0.2) .. image-sg:: /auto_calibration/least_squares_and_gaussian_calibration/images/sphx_glr_plot_calibration_flooding_005.png :alt: plot calibration flooding :srcset: /auto_calibration/least_squares_and_gaussian_calibration/images/sphx_glr_plot_calibration_flooding_005.png :class: sphx-glr-single-img .. GENERATED FROM PYTHON SOURCE LINES 430-466 .. code-block:: Python def plotDistributionGridPDF(distribution): """ Plot the marginal and bi-dimensional iso-PDF on a grid. Parameters ---------- distribution : ot.Distribution The distribution. Returns ------- grid : ot.GridLayout(dimension, dimension) The grid of plots. """ dimension = distribution.getDimension() grid = ot.GridLayout(dimension, dimension) for i in range(dimension): for j in range(dimension): if i == j: distributionI = distribution.getMarginal([i]) graph = distributionI.drawPDF() else: distributionJI = distribution.getMarginal([j, i]) graph = distributionJI.drawPDF() graph.setLegends([""]) graph.setTitle("") if i < dimension - 1: graph.setXTitle("") if j > 0: graph.setYTitle("") grid.setGraph(i, j, graph) grid.setTitle("Iso-PDF values") return grid .. GENERATED FROM PYTHON SOURCE LINES 467-468 Plot the PDF values of the distribution of the optimum parameters. .. GENERATED FROM PYTHON SOURCE LINES 468-477 .. code-block:: Python grid = plotDistributionGridPDF(thetaPosterior) view = otv.View( grid, figure_kw={"figsize": (6.0, 6.0)}, legend_kw={"bbox_to_anchor": (1.0, 1.0), "loc": "upper left"}, ) plot_space = 0.5 plt.subplots_adjust(wspace=plot_space, hspace=plot_space) .. image-sg:: /auto_calibration/least_squares_and_gaussian_calibration/images/sphx_glr_plot_calibration_flooding_006.png :alt: Iso-PDF values :srcset: /auto_calibration/least_squares_and_gaussian_calibration/images/sphx_glr_plot_calibration_flooding_006.png :class: sphx-glr-single-img .. GENERATED FROM PYTHON SOURCE LINES 478-480 Gaussian linear calibration --------------------------- .. GENERATED FROM PYTHON SOURCE LINES 482-483 The standard deviation of the observations. .. GENERATED FROM PYTHON SOURCE LINES 483-485 .. code-block:: Python sigmaH = 0.5 # (m^2) .. GENERATED FROM PYTHON SOURCE LINES 486-487 Define the covariance matrix of the output Y of the model. .. GENERATED FROM PYTHON SOURCE LINES 487-490 .. code-block:: Python errorCovariance = ot.CovarianceMatrix(1) errorCovariance[0, 0] = sigmaH**2 .. GENERATED FROM PYTHON SOURCE LINES 491-492 Define the covariance matrix of the parameters :math:`\theta` to calibrate. .. GENERATED FROM PYTHON SOURCE LINES 492-502 .. code-block:: Python sigmaKs = 5.0 sigmaZv = 1.0 sigmaZm = 1.0 # sigma = ot.CovarianceMatrix(3) sigma[0, 0] = sigmaKs**2 sigma[1, 1] = sigmaZv**2 sigma[2, 2] = sigmaZm**2 print(sigma) .. rst-class:: sphx-glr-script-out .. code-block:: none [[ 25 0 0 ] [ 0 1 0 ] [ 0 0 1 ]] .. GENERATED FROM PYTHON SOURCE LINES 503-505 The :class:`~openturns.GaussianLinearCalibration` class performs Gaussian linear calibration by linearizing the model in the neighbourhood of the prior. .. GENERATED FROM PYTHON SOURCE LINES 505-509 .. code-block:: Python algo = ot.GaussianLinearCalibration( mycf, Qobs, Hobs, thetaPrior, sigma, errorCovariance, "SVD" ) .. GENERATED FROM PYTHON SOURCE LINES 510-512 The :meth:`~openturns.GaussianLinearCalibration.run` method computes the solution of the problem. .. GENERATED FROM PYTHON SOURCE LINES 512-515 .. code-block:: Python algo.run() calibrationResult = algo.getResult() .. GENERATED FROM PYTHON SOURCE LINES 516-520 Analysis of the results ----------------------- The :meth:`~openturns.CalibrationResult.getParameterMAP` method returns the maximum of the posterior distribution of :math:`\theta`. .. GENERATED FROM PYTHON SOURCE LINES 520-523 .. code-block:: Python thetaMAP = calibrationResult.getParameterMAP() print(thetaMAP) .. rst-class:: sphx-glr-script-out .. code-block:: none [24.4058,48.1188,51.8812] .. GENERATED FROM PYTHON SOURCE LINES 524-528 .. code-block:: Python graph = calibrationResult.drawObservationsVsInputs() graph.setLegendPosition("upper left") view = otv.View(graph) .. image-sg:: /auto_calibration/least_squares_and_gaussian_calibration/images/sphx_glr_plot_calibration_flooding_007.png :alt: plot calibration flooding :srcset: /auto_calibration/least_squares_and_gaussian_calibration/images/sphx_glr_plot_calibration_flooding_007.png :class: sphx-glr-single-img .. GENERATED FROM PYTHON SOURCE LINES 529-531 We see that the output of the model after calibration is closer to the observations. .. GENERATED FROM PYTHON SOURCE LINES 533-536 .. code-block:: Python graph = calibrationResult.drawObservationsVsPredictions() view = otv.View(graph) .. image-sg:: /auto_calibration/least_squares_and_gaussian_calibration/images/sphx_glr_plot_calibration_flooding_008.png :alt: plot calibration flooding :srcset: /auto_calibration/least_squares_and_gaussian_calibration/images/sphx_glr_plot_calibration_flooding_008.png :class: sphx-glr-single-img .. GENERATED FROM PYTHON SOURCE LINES 537-538 In this case, the fit is satisfactory after calibration. .. GENERATED FROM PYTHON SOURCE LINES 540-544 .. code-block:: Python graph = calibrationResult.drawResiduals() graph.setLegendPosition("upper left") view = otv.View(graph) .. image-sg:: /auto_calibration/least_squares_and_gaussian_calibration/images/sphx_glr_plot_calibration_flooding_009.png :alt: Residual analysis :srcset: /auto_calibration/least_squares_and_gaussian_calibration/images/sphx_glr_plot_calibration_flooding_009.png :class: sphx-glr-single-img .. GENERATED FROM PYTHON SOURCE LINES 545-547 We see that the distribution of the residual is centered on zero after calibration. .. GENERATED FROM PYTHON SOURCE LINES 549-551 The :meth:`~openturns.CalibrationResult.getParameterPosterior` method returns the posterior normal distribution of :math:`\theta`. .. GENERATED FROM PYTHON SOURCE LINES 553-556 .. code-block:: Python distributionPosterior = calibrationResult.getParameterPosterior() print(distributionPosterior) .. rst-class:: sphx-glr-script-out .. code-block:: none Normal(mu = [24.4058,48.1188,51.8812], sigma = [4.09428,0.818856,0.818856], R = [[ 1 0.491369 -0.491369 ] [ 0.491369 1 0.491369 ] [ -0.491369 0.491369 1 ]]) .. GENERATED FROM PYTHON SOURCE LINES 557-558 We can compute a 95% credibility interval of the parameter :math:`\theta^\star`. .. GENERATED FROM PYTHON SOURCE LINES 558-564 .. code-block:: Python print( distributionPosterior.computeBilateralConfidenceIntervalWithMarginalProbability( 0.95 )[0] ) .. rst-class:: sphx-glr-script-out .. code-block:: none [14.8044, 34.0072] [46.1986, 50.0391] [49.9609, 53.8014] .. GENERATED FROM PYTHON SOURCE LINES 565-567 We see that there is a large uncertainty on the value of the parameter :math:`K_s` which can be as small as :math:`14` and as large as :math:`34`. .. GENERATED FROM PYTHON SOURCE LINES 569-570 We can compare the prior and posterior distributions of the marginals of :math:`\theta`. .. GENERATED FROM PYTHON SOURCE LINES 572-580 .. code-block:: Python graph = calibrationResult.drawParameterDistributions() view = otv.View( graph, figure_kw={"figsize": (8.0, 3.0)}, legend_kw={"bbox_to_anchor": (1.0, 1.0), "loc": "upper left"}, ) plt.subplots_adjust(right=0.8, bottom=0.2) .. image-sg:: /auto_calibration/least_squares_and_gaussian_calibration/images/sphx_glr_plot_calibration_flooding_010.png :alt: plot calibration flooding :srcset: /auto_calibration/least_squares_and_gaussian_calibration/images/sphx_glr_plot_calibration_flooding_010.png :class: sphx-glr-single-img .. GENERATED FROM PYTHON SOURCE LINES 581-588 The two distributions are different, which shows that the calibration is sensitive to the observations (if the observations were not sensitive, the two distributions were superimposed). Moreover, the two distributions are quite close, which implies that the prior distribution has played a role in the calibration (otherwise the two distributions would be completely different, indicating that only the observations were taken into account). .. GENERATED FROM PYTHON SOURCE LINES 590-591 Plot the PDF values of the distribution of the optimum parameters. .. GENERATED FROM PYTHON SOURCE LINES 591-600 .. code-block:: Python grid = plotDistributionGridPDF(thetaPosterior) view = otv.View( grid, figure_kw={"figsize": (6.0, 6.0)}, legend_kw={"bbox_to_anchor": (1.0, 1.0), "loc": "upper left"}, ) plot_space = 0.5 plt.subplots_adjust(wspace=plot_space, hspace=plot_space) .. image-sg:: /auto_calibration/least_squares_and_gaussian_calibration/images/sphx_glr_plot_calibration_flooding_011.png :alt: Iso-PDF values :srcset: /auto_calibration/least_squares_and_gaussian_calibration/images/sphx_glr_plot_calibration_flooding_011.png :class: sphx-glr-single-img .. GENERATED FROM PYTHON SOURCE LINES 601-605 Gaussian nonlinear calibration ------------------------------ The :class:`~openturns.GaussianNonLinearCalibration` class performs Gaussian nonlinear calibration. .. GENERATED FROM PYTHON SOURCE LINES 605-609 .. code-block:: Python algo = ot.GaussianNonLinearCalibration( mycf, Qobs, Hobs, thetaPrior, sigma, errorCovariance ) .. GENERATED FROM PYTHON SOURCE LINES 610-612 The :meth:`~openturns.GaussianNonLinearCalibration.run` method computes the solution of the problem. .. GENERATED FROM PYTHON SOURCE LINES 612-615 .. code-block:: Python algo.run() calibrationResult = algo.getResult() .. GENERATED FROM PYTHON SOURCE LINES 616-620 Analysis of the results ----------------------- The :meth:`~openturns.CalibrationResult.getParameterMAP` method returns the maximum of the posterior distribution of :math:`\theta`. .. GENERATED FROM PYTHON SOURCE LINES 620-623 .. code-block:: Python thetaMAP = calibrationResult.getParameterMAP() print(thetaMAP) .. rst-class:: sphx-glr-script-out .. code-block:: none [29.4528,47.7599,52.2401] .. GENERATED FROM PYTHON SOURCE LINES 624-628 .. code-block:: Python graph = calibrationResult.drawObservationsVsInputs() graph.setLegendPosition("upper left") view = otv.View(graph) .. image-sg:: /auto_calibration/least_squares_and_gaussian_calibration/images/sphx_glr_plot_calibration_flooding_012.png :alt: plot calibration flooding :srcset: /auto_calibration/least_squares_and_gaussian_calibration/images/sphx_glr_plot_calibration_flooding_012.png :class: sphx-glr-single-img .. GENERATED FROM PYTHON SOURCE LINES 629-631 We see that the output of the model after calibration is in the middle of the observations: the calibration seems correct. .. GENERATED FROM PYTHON SOURCE LINES 633-636 .. code-block:: Python graph = calibrationResult.drawObservationsVsPredictions() view = otv.View(graph) .. image-sg:: /auto_calibration/least_squares_and_gaussian_calibration/images/sphx_glr_plot_calibration_flooding_013.png :alt: plot calibration flooding :srcset: /auto_calibration/least_squares_and_gaussian_calibration/images/sphx_glr_plot_calibration_flooding_013.png :class: sphx-glr-single-img .. GENERATED FROM PYTHON SOURCE LINES 637-638 The fit is excellent after calibration. Indeed, the cloud of points after calibration is on the diagonal. .. GENERATED FROM PYTHON SOURCE LINES 640-648 .. code-block:: Python graph = calibrationResult.drawResiduals() view = otv.View( graph, figure_kw={"figsize": (8.0, 4.0)}, legend_kw={"bbox_to_anchor": (1.0, 1.0), "loc": "upper left"}, ) plt.subplots_adjust(right=0.6) .. image-sg:: /auto_calibration/least_squares_and_gaussian_calibration/images/sphx_glr_plot_calibration_flooding_014.png :alt: Residual analysis :srcset: /auto_calibration/least_squares_and_gaussian_calibration/images/sphx_glr_plot_calibration_flooding_014.png :class: sphx-glr-single-img .. GENERATED FROM PYTHON SOURCE LINES 649-651 We see that the distribution of the residual is centered on zero. This is a proof that the calibration did perform correctly. .. GENERATED FROM PYTHON SOURCE LINES 653-655 The :meth:`~openturns.CalibrationResult.getParameterPosterior` method returns the posterior normal distribution of :math:`\theta`. .. GENERATED FROM PYTHON SOURCE LINES 655-657 .. code-block:: Python distributionPosterior = calibrationResult.getParameterPosterior() .. GENERATED FROM PYTHON SOURCE LINES 658-659 We can compute a 95% credibility interval of the parameter :math:`\theta^\star`. .. GENERATED FROM PYTHON SOURCE LINES 659-665 .. code-block:: Python print( distributionPosterior.computeBilateralConfidenceIntervalWithMarginalProbability( 0.95 )[0] ) .. rst-class:: sphx-glr-script-out .. code-block:: none [29.9114, 30.9625] [47.5676, 47.7] [52.3, 52.4324] .. GENERATED FROM PYTHON SOURCE LINES 666-667 We see that there is a small uncertainty on the value of all parameters. .. GENERATED FROM PYTHON SOURCE LINES 669-670 We can compare the prior and posterior distributions of the marginals of :math:`\theta`. .. GENERATED FROM PYTHON SOURCE LINES 670-678 .. code-block:: Python graph = calibrationResult.drawParameterDistributions() view = otv.View( graph, figure_kw={"figsize": (8.0, 3.0)}, legend_kw={"bbox_to_anchor": (1.0, 1.0), "loc": "upper left"}, ) plt.subplots_adjust(right=0.8, bottom=0.2) .. image-sg:: /auto_calibration/least_squares_and_gaussian_calibration/images/sphx_glr_plot_calibration_flooding_015.png :alt: plot calibration flooding :srcset: /auto_calibration/least_squares_and_gaussian_calibration/images/sphx_glr_plot_calibration_flooding_015.png :class: sphx-glr-single-img .. GENERATED FROM PYTHON SOURCE LINES 679-681 The two distributions are very different, with a spiky posterior distribution. This shows that the calibration is very sensitive to the observations. .. GENERATED FROM PYTHON SOURCE LINES 683-684 Plot the PDF values of the distribution of the optimum parameters. .. GENERATED FROM PYTHON SOURCE LINES 684-693 .. code-block:: Python grid = plotDistributionGridPDF(thetaPosterior) view = otv.View( grid, figure_kw={"figsize": (6.0, 6.0)}, legend_kw={"bbox_to_anchor": (1.0, 1.0), "loc": "upper left"}, ) plot_space = 0.5 plt.subplots_adjust(wspace=plot_space, hspace=plot_space) .. image-sg:: /auto_calibration/least_squares_and_gaussian_calibration/images/sphx_glr_plot_calibration_flooding_016.png :alt: Iso-PDF values :srcset: /auto_calibration/least_squares_and_gaussian_calibration/images/sphx_glr_plot_calibration_flooding_016.png :class: sphx-glr-single-img .. GENERATED FROM PYTHON SOURCE LINES 694-708 Tuning the posterior distribution estimation -------------------------------------------- The "GaussianNonLinearCalibration-BootstrapSize" key of the :class:`~openturns.ResourceMap` controls the posterior distribution estimation. * If "GaussianNonLinearCalibration-BootstrapSize" > 0 (by default it is equal to 100), then a bootstrap resample algorithm is used to see the dispersion of the MAP estimator. This allows one to see the variability of the estimator with respect to the finite noisy observation sample. * If "GaussianNonLinearCalibration-BootstrapSize" is zero, then the Gaussian linear calibration estimator is used (i.e. the :class:`~openturns.GaussianLinearCalibration` class) at the optimum. This is called the Laplace approximation. .. GENERATED FROM PYTHON SOURCE LINES 710-713 The default value of the key is nonzero, meaning that bootstrap is used. This can be costly in some cases, because it requires to repeat the optimization several times. .. GENERATED FROM PYTHON SOURCE LINES 713-715 .. code-block:: Python print(ot.ResourceMap.GetAsUnsignedInteger("GaussianNonLinearCalibration-BootstrapSize")) .. rst-class:: sphx-glr-script-out .. code-block:: none 100 .. GENERATED FROM PYTHON SOURCE LINES 716-718 We must configure the key before creating the object (otherwise changing the parameter does not change the result). .. GENERATED FROM PYTHON SOURCE LINES 718-732 .. code-block:: Python ot.ResourceMap.SetAsUnsignedInteger("GaussianNonLinearCalibration-BootstrapSize", 0) algo = ot.GaussianNonLinearCalibration( mycf, Qobs, Hobs, thetaPrior, sigma, errorCovariance ) algo.run() calibrationResult = algo.getResult() graph = calibrationResult.drawParameterDistributions() view = otv.View( graph, figure_kw={"figsize": (8.0, 3.0)}, legend_kw={"bbox_to_anchor": (1.0, 1.0), "loc": "upper left"}, ) plt.subplots_adjust(right=0.8, bottom=0.2) .. image-sg:: /auto_calibration/least_squares_and_gaussian_calibration/images/sphx_glr_plot_calibration_flooding_017.png :alt: plot calibration flooding :srcset: /auto_calibration/least_squares_and_gaussian_calibration/images/sphx_glr_plot_calibration_flooding_017.png :class: sphx-glr-single-img .. GENERATED FROM PYTHON SOURCE LINES 733-734 As we can see, this does not change much the posterior distribution, which remains spiky. .. GENERATED FROM PYTHON SOURCE LINES 736-737 Plot the PDF values of the distribution of the optimum parameters. .. GENERATED FROM PYTHON SOURCE LINES 737-748 .. code-block:: Python grid = plotDistributionGridPDF(thetaPosterior) view = otv.View( grid, figure_kw={"figsize": (6.0, 6.0)}, legend_kw={"bbox_to_anchor": (1.0, 1.0), "loc": "upper left"}, ) plot_space = 0.5 plt.subplots_adjust(wspace=plot_space, hspace=plot_space) otv.View.ShowAll() .. image-sg:: /auto_calibration/least_squares_and_gaussian_calibration/images/sphx_glr_plot_calibration_flooding_018.png :alt: Iso-PDF values :srcset: /auto_calibration/least_squares_and_gaussian_calibration/images/sphx_glr_plot_calibration_flooding_018.png :class: sphx-glr-single-img .. GENERATED FROM PYTHON SOURCE LINES 749-750 Reset default settings .. GENERATED FROM PYTHON SOURCE LINES 750-751 .. code-block:: Python ot.ResourceMap.Reload() .. rst-class:: sphx-glr-timing **Total running time of the script:** (0 minutes 8.856 seconds) .. _sphx_glr_download_auto_calibration_least_squares_and_gaussian_calibration_plot_calibration_flooding.py: .. only:: html .. container:: sphx-glr-footer sphx-glr-footer-example .. container:: sphx-glr-download sphx-glr-download-jupyter :download:`Download Jupyter notebook: plot_calibration_flooding.ipynb ` .. container:: sphx-glr-download sphx-glr-download-python :download:`Download Python source code: plot_calibration_flooding.py `