.. 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 <sphx_glr_download_auto_calibration_least_squares_and_gaussian_calibration_plot_calibration_flooding.py>`
        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 <use-case-flood-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
</auto_calibration/least_squares_and_gaussian_calibration/plot_calibration_chaboche>` 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
</auto_calibration/least_squares_and_gaussian_calibration/plot_generate_flooding>`.

.. GENERATED FROM PYTHON SOURCE LINES 21-49

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<use-case-flood-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 49-58

.. code-block:: Python


    from openturns.usecases import flood_model
    from matplotlib import pyplot as plt
    import openturns.viewer as otv
    import numpy as np
    import openturns as ot

    ot.ResourceMap.SetAsUnsignedInteger("Normal-SmallDimension", 1)








.. 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

    <div class="output_subarea output_html rendered_html output_result">
    <table>
      <tr><td></td><th>Q ($m^3/s$)</th><th>H (m)</th></tr>
      <tr><th>0</th><td>130</td><td>0.59</td></tr>
      <tr><th>1</th><td>530</td><td>1.33</td></tr>
      <tr><th>2</th><td>960</td><td>2.03</td></tr>
      <tr><th>3</th><td>1400</td><td>2.72</td></tr>
      <tr><th>4</th><td>1830</td><td>2.83</td></tr>
      <tr><th>5</th><td>2260</td><td>3.5</td></tr>
      <tr><th>6</th><td>2700</td><td>3.82</td></tr>
      <tr><th>7</th><td>3130</td><td>4.36</td></tr>
      <tr><th>8</th><td>3560</td><td>4.63</td></tr>
      <tr><th>9</th><td>4010</td><td>4.96</td></tr>
    </table>
    </div>
    <br />
    <br />

.. 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        |
+-------+-----------+
| 2     | 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.svg
   :alt: Observations
   :srcset: /auto_calibration/least_squares_and_gaussian_calibration/images/sphx_glr_plot_calibration_flooding_001.svg
   :class: sphx-glr-single-img





.. GENERATED FROM PYTHON SOURCE LINES 205-208

We 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,-8.8645e+23,-8.8645e+23]
    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,-8.8645e+23,-8.8645e+23], sigma = [3.2656e+27,2.30913e+33,2.30913e+33], R = [[  1            5.71128e-25 -5.71128e-25 ]
     [  5.71128e-25  1            1           ]
     [ -5.71128e-25  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

    [-7.2902e+27, 7.2902e+27]
    [-5.15495e+33, 5.15495e+33]
    [-5.15495e+33, 5.15495e+33]




.. 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,1.08063e-25]




.. 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.2562, 27.922]
    [46.9414, 47.2226]
    [52.7774, 53.0586]




.. 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.svg
   :alt: plot calibration flooding
   :srcset: /auto_calibration/least_squares_and_gaussian_calibration/images/sphx_glr_plot_calibration_flooding_002.svg
   :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.svg
   :alt: plot calibration flooding
   :srcset: /auto_calibration/least_squares_and_gaussian_calibration/images/sphx_glr_plot_calibration_flooding_003.svg
   :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.svg
   :alt: Residual analysis
   :srcset: /auto_calibration/least_squares_and_gaussian_calibration/images/sphx_glr_plot_calibration_flooding_004.svg
   :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.svg
   :alt: plot calibration flooding
   :srcset: /auto_calibration/least_squares_and_gaussian_calibration/images/sphx_glr_plot_calibration_flooding_005.svg
   :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-476

.. code-block:: Python

    grid = plotDistributionGridPDF(thetaPosterior)
    view = otv.View(
        grid,
        figure_kw={"figsize": (9.0, 6.0)},
        legend_kw={"bbox_to_anchor": (1.0, 1.0), "loc": "upper left"},
    )
    plt.subplots_adjust(wspace=1.0, hspace=0.5)




.. image-sg:: /auto_calibration/least_squares_and_gaussian_calibration/images/sphx_glr_plot_calibration_flooding_006.svg
   :alt: Iso-PDF values
   :srcset: /auto_calibration/least_squares_and_gaussian_calibration/images/sphx_glr_plot_calibration_flooding_006.svg
   :class: sphx-glr-single-img





.. GENERATED FROM PYTHON SOURCE LINES 477-479

Gaussian linear calibration
---------------------------

.. GENERATED FROM PYTHON SOURCE LINES 481-482

The standard deviation of the observations.

.. GENERATED FROM PYTHON SOURCE LINES 482-484

.. code-block:: Python

    sigmaH = 0.5  # (m^2)








.. GENERATED FROM PYTHON SOURCE LINES 485-486

Define the covariance matrix of the output Y of the model.

.. GENERATED FROM PYTHON SOURCE LINES 486-489

.. code-block:: Python

    errorCovariance = ot.CovarianceMatrix(1)
    errorCovariance[0, 0] = sigmaH**2








.. GENERATED FROM PYTHON SOURCE LINES 490-491

Define the covariance matrix of the parameters :math:`\theta` to calibrate.

.. GENERATED FROM PYTHON SOURCE LINES 491-501

.. 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 502-504

The :class:`~openturns.GaussianLinearCalibration` class performs Gaussian
linear calibration by linearizing the model in the neighbourhood of the prior.

.. GENERATED FROM PYTHON SOURCE LINES 504-508

.. code-block:: Python

    algo = ot.GaussianLinearCalibration(
        mycf, Qobs, Hobs, thetaPrior, sigma, errorCovariance, "SVD"
    )








.. GENERATED FROM PYTHON SOURCE LINES 509-511

The :meth:`~openturns.GaussianLinearCalibration.run` method computes
the solution of the problem.

.. GENERATED FROM PYTHON SOURCE LINES 511-514

.. code-block:: Python

    algo.run()
    calibrationResult = algo.getResult()








.. GENERATED FROM PYTHON SOURCE LINES 515-519

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 519-522

.. 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 523-527

.. 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.svg
   :alt: plot calibration flooding
   :srcset: /auto_calibration/least_squares_and_gaussian_calibration/images/sphx_glr_plot_calibration_flooding_007.svg
   :class: sphx-glr-single-img





.. GENERATED FROM PYTHON SOURCE LINES 528-530

We see that the output of the model after calibration is closer to the
observations.

.. GENERATED FROM PYTHON SOURCE LINES 532-535

.. 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.svg
   :alt: plot calibration flooding
   :srcset: /auto_calibration/least_squares_and_gaussian_calibration/images/sphx_glr_plot_calibration_flooding_008.svg
   :class: sphx-glr-single-img





.. GENERATED FROM PYTHON SOURCE LINES 536-537

In this case, the fit is satisfactory after calibration.

.. GENERATED FROM PYTHON SOURCE LINES 539-543

.. 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.svg
   :alt: Residual analysis
   :srcset: /auto_calibration/least_squares_and_gaussian_calibration/images/sphx_glr_plot_calibration_flooding_009.svg
   :class: sphx-glr-single-img





.. GENERATED FROM PYTHON SOURCE LINES 544-546

We see that the distribution of the residual is centered on zero after
calibration.

.. GENERATED FROM PYTHON SOURCE LINES 548-550

The :meth:`~openturns.CalibrationResult.getParameterPosterior` method
returns the posterior normal distribution of :math:`\theta`.

.. GENERATED FROM PYTHON SOURCE LINES 552-555

.. 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 556-557

We can compute a 95% credibility interval of the parameter :math:`\theta^\star`.

.. GENERATED FROM PYTHON SOURCE LINES 557-563

.. code-block:: Python

    print(
        distributionPosterior.computeBilateralConfidenceIntervalWithMarginalProbability(
            0.95
        )[0]
    )





.. rst-class:: sphx-glr-script-out

 .. code-block:: none

    [14.8009, 34.0107]
    [46.1979, 50.0398]
    [49.9602, 53.8021]




.. GENERATED FROM PYTHON SOURCE LINES 564-566

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 568-569

We can compare the prior and posterior distributions of the marginals of :math:`\theta`.

.. GENERATED FROM PYTHON SOURCE LINES 571-579

.. 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.svg
   :alt: plot calibration flooding
   :srcset: /auto_calibration/least_squares_and_gaussian_calibration/images/sphx_glr_plot_calibration_flooding_010.svg
   :class: sphx-glr-single-img





.. GENERATED FROM PYTHON SOURCE LINES 580-587

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 589-590

Plot the PDF values of the distribution of the optimum parameters.

.. GENERATED FROM PYTHON SOURCE LINES 590-598

.. code-block:: Python

    grid = plotDistributionGridPDF(thetaPosterior)
    view = otv.View(
        grid,
        figure_kw={"figsize": (9.0, 6.0)},
        legend_kw={"bbox_to_anchor": (1.0, 1.0), "loc": "upper left"},
    )
    plt.subplots_adjust(wspace=1.0, hspace=0.5)




.. image-sg:: /auto_calibration/least_squares_and_gaussian_calibration/images/sphx_glr_plot_calibration_flooding_011.svg
   :alt: Iso-PDF values
   :srcset: /auto_calibration/least_squares_and_gaussian_calibration/images/sphx_glr_plot_calibration_flooding_011.svg
   :class: sphx-glr-single-img





.. GENERATED FROM PYTHON SOURCE LINES 599-603

Gaussian nonlinear calibration
------------------------------
The :class:`~openturns.GaussianNonLinearCalibration` class performs Gaussian
nonlinear calibration.

.. GENERATED FROM PYTHON SOURCE LINES 603-607

.. code-block:: Python

    algo = ot.GaussianNonLinearCalibration(
        mycf, Qobs, Hobs, thetaPrior, sigma, errorCovariance
    )








.. GENERATED FROM PYTHON SOURCE LINES 608-610

The :meth:`~openturns.GaussianNonLinearCalibration.run` method computes the
solution of the problem.

.. GENERATED FROM PYTHON SOURCE LINES 610-613

.. code-block:: Python

    algo.run()
    calibrationResult = algo.getResult()








.. GENERATED FROM PYTHON SOURCE LINES 614-618

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 618-621

.. 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 622-626

.. 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.svg
   :alt: plot calibration flooding
   :srcset: /auto_calibration/least_squares_and_gaussian_calibration/images/sphx_glr_plot_calibration_flooding_012.svg
   :class: sphx-glr-single-img





.. GENERATED FROM PYTHON SOURCE LINES 627-629

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 631-634

.. 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.svg
   :alt: plot calibration flooding
   :srcset: /auto_calibration/least_squares_and_gaussian_calibration/images/sphx_glr_plot_calibration_flooding_013.svg
   :class: sphx-glr-single-img





.. GENERATED FROM PYTHON SOURCE LINES 635-636

The fit is excellent after calibration. Indeed, the cloud of points after calibration is on the diagonal.

.. GENERATED FROM PYTHON SOURCE LINES 638-646

.. 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.svg
   :alt: Residual analysis
   :srcset: /auto_calibration/least_squares_and_gaussian_calibration/images/sphx_glr_plot_calibration_flooding_014.svg
   :class: sphx-glr-single-img





.. GENERATED FROM PYTHON SOURCE LINES 647-649

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 651-653

The :meth:`~openturns.CalibrationResult.getParameterPosterior` method
returns the posterior normal distribution of :math:`\theta`.

.. GENERATED FROM PYTHON SOURCE LINES 653-655

.. code-block:: Python

    distributionPosterior = calibrationResult.getParameterPosterior()








.. GENERATED FROM PYTHON SOURCE LINES 656-657

We can compute a 95% credibility interval of the parameter :math:`\theta^\star`.

.. GENERATED FROM PYTHON SOURCE LINES 657-663

.. code-block:: Python

    print(
        distributionPosterior.computeBilateralConfidenceIntervalWithMarginalProbability(
            0.95
        )[0]
    )





.. rst-class:: sphx-glr-script-out

 .. code-block:: none

    [29.9146, 30.8866]
    [47.5769, 47.6996]
    [52.3004, 52.4231]




.. GENERATED FROM PYTHON SOURCE LINES 664-665

We see that there is a small uncertainty on the value of all parameters.

.. GENERATED FROM PYTHON SOURCE LINES 667-668

We can compare the prior and posterior distributions of the marginals of :math:`\theta`.

.. GENERATED FROM PYTHON SOURCE LINES 668-676

.. 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.svg
   :alt: plot calibration flooding
   :srcset: /auto_calibration/least_squares_and_gaussian_calibration/images/sphx_glr_plot_calibration_flooding_015.svg
   :class: sphx-glr-single-img





.. GENERATED FROM PYTHON SOURCE LINES 677-679

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 681-682

Plot the PDF values of the distribution of the optimum parameters.

.. GENERATED FROM PYTHON SOURCE LINES 682-690

.. code-block:: Python

    grid = plotDistributionGridPDF(thetaPosterior)
    view = otv.View(
        grid,
        figure_kw={"figsize": (9.0, 6.0)},
        legend_kw={"bbox_to_anchor": (1.0, 1.0), "loc": "upper left"},
    )
    plt.subplots_adjust(wspace=1.0, hspace=0.5)




.. image-sg:: /auto_calibration/least_squares_and_gaussian_calibration/images/sphx_glr_plot_calibration_flooding_016.svg
   :alt: Iso-PDF values
   :srcset: /auto_calibration/least_squares_and_gaussian_calibration/images/sphx_glr_plot_calibration_flooding_016.svg
   :class: sphx-glr-single-img





.. GENERATED FROM PYTHON SOURCE LINES 691-705

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 707-710

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 710-712

.. code-block:: Python

    print(ot.ResourceMap.GetAsUnsignedInteger("GaussianNonLinearCalibration-BootstrapSize"))





.. rst-class:: sphx-glr-script-out

 .. code-block:: none

    100




.. GENERATED FROM PYTHON SOURCE LINES 713-715

We must configure the key before creating the object (otherwise changing
the parameter does not change the result).

.. GENERATED FROM PYTHON SOURCE LINES 715-729

.. 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.svg
   :alt: plot calibration flooding
   :srcset: /auto_calibration/least_squares_and_gaussian_calibration/images/sphx_glr_plot_calibration_flooding_017.svg
   :class: sphx-glr-single-img





.. GENERATED FROM PYTHON SOURCE LINES 730-731

As we can see, this does not change much the posterior distribution, which remains spiky.

.. GENERATED FROM PYTHON SOURCE LINES 733-734

Plot the PDF values of the distribution of the optimum parameters.

.. GENERATED FROM PYTHON SOURCE LINES 734-742

.. code-block:: Python

    grid = plotDistributionGridPDF(thetaPosterior)
    view = otv.View(
        grid,
        figure_kw={"figsize": (9.0, 6.0)},
        legend_kw={"bbox_to_anchor": (1.0, 1.0), "loc": "upper left"},
    )
    plt.subplots_adjust(wspace=1.0, hspace=0.5)




.. image-sg:: /auto_calibration/least_squares_and_gaussian_calibration/images/sphx_glr_plot_calibration_flooding_018.svg
   :alt: Iso-PDF values
   :srcset: /auto_calibration/least_squares_and_gaussian_calibration/images/sphx_glr_plot_calibration_flooding_018.svg
   :class: sphx-glr-single-img





.. GENERATED FROM PYTHON SOURCE LINES 743-744

.. code-block:: Python

    otv.View.ShowAll()








.. _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 <plot_calibration_flooding.ipynb>`

    .. container:: sphx-glr-download sphx-glr-download-python

      :download:`Download Python source code: plot_calibration_flooding.py <plot_calibration_flooding.py>`

    .. container:: sphx-glr-download sphx-glr-download-zip

      :download:`Download zipped: plot_calibration_flooding.zip <plot_calibration_flooding.zip>`