.. 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_quickstart.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_quickstart.py: Calibrate a parametric model: a quick-start guide to calibration ================================================================ In this example we present the calibration of a parametric module. To do this, we show how to define the observations. Then we create a Python function and create the parametric model that is to be calibrated. Finally, we calibrate the parameters of the model using least squares and we validate the hypotheses of the method. Please read :ref:`code_calibration` for more details on code calibration and least squares. In this example we are interested in the calibration of the :ref:`flooding model `. Once the reader has mastered this example, the :doc:`calibration of the Chaboche mechanical model ` may be considered to make an in-depth study of these algorithms. The goal of this script is to be relatively small, so please consider reading the other examples if this is relevant. .. GENERATED FROM PYTHON SOURCE LINES 21-44 Parameters to calibrate and observations ---------------------------------------- The parametric model is: .. math:: H = a Q^b. The vector of parameters to calibrate is: .. math:: \theta = (a,b). 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. We choose to calibrate this model using non linear least squares, because this is a relatively flexible method. .. GENERATED FROM PYTHON SOURCE LINES 46-53 .. code-block:: Python import openturns as ot import openturns.viewer as otv from openturns.usecases import flood_model ot.Log.Show(ot.Log.NONE) .. GENERATED FROM PYTHON SOURCE LINES 54-62 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 64-71 .. 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 72-81 Create the physical model ------------------------- We define the model :math:`g` which has 3 inputs and one output `H`. This model has a parametric shape that may correspond to the data, which has some power shape. In the model, the parameters are `a` and `b`, the input is `Q` and the output is `H`: - `Q` : the flowrate of the river, - `a`, `b` : the parameters. .. GENERATED FROM PYTHON SOURCE LINES 84-96 .. code-block:: Python def functionSimpleFlooding(X): Q, a, b = X H = a * Q**b return [H] g = ot.PythonFunction(3, 1, functionSimpleFlooding) g = ot.MemoizeFunction(g) g.setInputDescription(["Q", "a", "b"]) g.setOutputDescription(["H"]) .. GENERATED FROM PYTHON SOURCE LINES 97-99 Setting the calibration parameters ---------------------------------- .. GENERATED FROM PYTHON SOURCE LINES 101-104 Define the value of the reference values of the :math:`\theta` parameter. There is no particular method to set these values: we used trial-and-error to see the order of magnitude of the parameters. .. GENERATED FROM PYTHON SOURCE LINES 106-110 .. code-block:: Python aInitial = 0.1 bInitial = 0.5 thetaPrior = [aInitial, bInitial] .. GENERATED FROM PYTHON SOURCE LINES 111-135 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 | a | +-------+----------------+ | 2 | b | +-------+----------------+ +-------+-----------+ | Index | Parameter | +=======+===========+ | ∅ | ∅ | +-------+-----------+ **Table 1.** Indices and names of the inputs and parameters of the physical model. .. GENERATED FROM PYTHON SOURCE LINES 137-140 .. 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,a,b] Physical Model Parameters: [] .. GENERATED FROM PYTHON SOURCE LINES 141-165 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 `a` and `b` which have the indices 1 and and 2 in the physical model. Please read :doc:`Create a parametric function ` for more details on this topic. +-------+----------------+ | Index | Input variable | +=======+================+ | 0 | Q | +-------+----------------+ +-------+-----------+ | Index | Parameter | +=======+===========+ | 0 | a | +-------+-----------+ | 1 | b | +-------+-----------+ **Table 2.** Indices and names of the inputs and parameters of the parametric model. .. GENERATED FROM PYTHON SOURCE LINES 168-171 The following statement create the calibrated function from the model. The calibrated parameters :math:`a`, :math:`b` are at indices 1, 2 in the inputs arguments of the model. .. GENERATED FROM PYTHON SOURCE LINES 173-176 .. code-block:: Python calibratedIndices = [1, 2] mycf = ot.ParametricFunction(g, calibratedIndices, thetaPrior) .. GENERATED FROM PYTHON SOURCE LINES 177-178 Plot the Y observations versus the X observations. .. GENERATED FROM PYTHON SOURCE LINES 180-194 .. code-block:: Python title = "Before calibration : a = %.4f, b = %.4f" % (aInitial, bInitial) graph = ot.Graph(title, "Q", "H", True) # Plot the model before calibration curve = mycf.draw(100.0, 4000.0).getDrawable(0) curve.setLegend("Model, before calibration") graph.add(curve) # Plot the noisy observations cloud = ot.Cloud(Qobs, Hobs) cloud.setLegend("Observations") 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_quickstart_001.png :alt: Before calibration : a = 0.1000, b = 0.5000 :srcset: /auto_calibration/least_squares_and_gaussian_calibration/images/sphx_glr_plot_calibration_quickstart_001.png :class: sphx-glr-single-img .. GENERATED FROM PYTHON SOURCE LINES 195-198 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 200-205 Calibration with non linear least squares ----------------------------------------- The :class:`~openturns.NonLinearLeastSquaresCalibration` class performs non linear least squares calibration by minimizing the squared Euclidian norm between the predictions and the observations. .. GENERATED FROM PYTHON SOURCE LINES 207-209 .. code-block:: Python algo = ot.NonLinearLeastSquaresCalibration(mycf, Qobs, Hobs, thetaPrior) .. GENERATED FROM PYTHON SOURCE LINES 210-212 The :meth:`~openturns.NonLinearLeastSquaresCalibration.run` method computes the solution of the problem. .. GENERATED FROM PYTHON SOURCE LINES 214-217 .. code-block:: Python algo.run() calibrationResult = algo.getResult() .. GENERATED FROM PYTHON SOURCE LINES 218-220 Analysis of the results ----------------------- .. GENERATED FROM PYTHON SOURCE LINES 222-224 The :meth:`~openturns.CalibrationResult.getParameterMAP` method returns the maximum of the posterior distribution of :math:`\theta`. .. GENERATED FROM PYTHON SOURCE LINES 226-230 .. code-block:: Python thetaMAP = calibrationResult.getParameterMAP() print(thetaMAP) .. rst-class:: sphx-glr-script-out .. code-block:: none [0.0258034,0.634465] .. GENERATED FROM PYTHON SOURCE LINES 231-233 In order to see if the parameters fit the data, we plot the input depending on the output before and after calibration. .. GENERATED FROM PYTHON SOURCE LINES 233-243 .. code-block:: Python # sphinx_gallery_thumbnail_number = 2 graph = calibrationResult.drawObservationsVsInputs() aEstimated, bEstimated = thetaMAP title = "After calibration : a = %.4f, b = %.4f" % (aEstimated, bEstimated) graph.setTitle(title) graph.setLegendPosition("upper left") view = otv.View(graph) .. image-sg:: /auto_calibration/least_squares_and_gaussian_calibration/images/sphx_glr_plot_calibration_quickstart_002.png :alt: After calibration : a = 0.0258, b = 0.6345 :srcset: /auto_calibration/least_squares_and_gaussian_calibration/images/sphx_glr_plot_calibration_quickstart_002.png :class: sphx-glr-single-img .. GENERATED FROM PYTHON SOURCE LINES 244-247 One of the hypotheses of the least squares method is that the residuals follow a normal distribution: the next cell checks if this hypothesis is satisfied here. .. GENERATED FROM PYTHON SOURCE LINES 249-253 .. 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_quickstart_003.png :alt: Residual analysis :srcset: /auto_calibration/least_squares_and_gaussian_calibration/images/sphx_glr_plot_calibration_quickstart_003.png :class: sphx-glr-single-img .. GENERATED FROM PYTHON SOURCE LINES 254-258 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. .. GENERATED FROM PYTHON SOURCE LINES 260-271 *Remark.* The logarithm of the height is: .. math:: \log(H) = \log(a) + b \log(Q). Hence, transforming the data on a logarithmic scale leads to a parametric model that is linear in the parameters. The parameters of this transformed model can be estimated using linear linear squares, which may lead to a significant improvement in terms of number of function evaluations. .. GENERATED FROM PYTHON SOURCE LINES 271-273 .. code-block:: Python otv.View.ShowAll() .. _sphx_glr_download_auto_calibration_least_squares_and_gaussian_calibration_plot_calibration_quickstart.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_quickstart.ipynb ` .. container:: sphx-glr-download sphx-glr-download-python :download:`Download Python source code: plot_calibration_quickstart.py `