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

    <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        |
+-------+-----------+
| 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-205

.. 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.setColors(ot.Drawable.BuildDefaultPalette(2))
    graph.setLegendPosition("topleft")
    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 206-209

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 211-216

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 216-218

.. code-block:: Python

    algo = ot.LinearLeastSquaresCalibration(mycf, Qobs, Hobs, thetaPrior, "SVD")








.. GENERATED FROM PYTHON SOURCE LINES 219-221

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

.. GENERATED FROM PYTHON SOURCE LINES 221-224

.. code-block:: Python

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








.. GENERATED FROM PYTHON SOURCE LINES 225-227

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

.. GENERATED FROM PYTHON SOURCE LINES 227-230

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

Print the true values of the parameters.

.. GENERATED FROM PYTHON SOURCE LINES 232-237

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

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 244-246

Diagnostic of the identification issue
--------------------------------------

.. GENERATED FROM PYTHON SOURCE LINES 248-252

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 254-257

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

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 263-269

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

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 275-282

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

The rank of the problem can be seen from the singular values of the Jacobian
matrix.

.. GENERATED FROM PYTHON SOURCE LINES 287-289

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

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 309-311

Conclusion of the linear least squares calibration
--------------------------------------------------

.. GENERATED FROM PYTHON SOURCE LINES 313-325

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 327-332

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 334-336

.. code-block:: Python

    algo = ot.NonLinearLeastSquaresCalibration(mycf, Qobs, Hobs, thetaPrior)








.. GENERATED FROM PYTHON SOURCE LINES 337-339

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

.. GENERATED FROM PYTHON SOURCE LINES 339-342

.. code-block:: Python

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








.. GENERATED FROM PYTHON SOURCE LINES 343-345

Analysis of the results
-----------------------

.. GENERATED FROM PYTHON SOURCE LINES 347-349

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

.. GENERATED FROM PYTHON SOURCE LINES 351-354

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

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 362-368

.. code-block:: Python

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





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

 .. code-block:: none

    100




.. GENERATED FROM PYTHON SOURCE LINES 369-372

.. code-block:: Python

    thetaPosterior = calibrationResult.getParameterPosterior()
    print(thetaPosterior.computeBilateralConfidenceIntervalWithMarginalProbability(0.95)[0])





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

 .. code-block:: none

    [27.2377, 28.1072]
    [46.863, 47.2304]
    [52.7696, 53.137]




.. GENERATED FROM PYTHON SOURCE LINES 373-375

In this case, the values of the parameters are quite accurately
computed.

.. GENERATED FROM PYTHON SOURCE LINES 377-379

Increase the default number of points in the plots.
This produces smoother spiky distributions.

.. GENERATED FROM PYTHON SOURCE LINES 379-381

.. code-block:: Python

    ot.ResourceMap.SetAsUnsignedInteger("Distribution-DefaultPointNumber", 300)








.. GENERATED FROM PYTHON SOURCE LINES 382-386

.. code-block:: Python

    graph = calibrationResult.drawObservationsVsInputs()
    graph.setLegendPosition("topleft")
    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 387-389

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 391-394

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

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 399-402

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

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 405-411

.. code-block:: Python


    # sphinx_gallery_thumbnail_number = 5
    graph = calibrationResult.drawResiduals()
    graph.setLegendPosition("topleft")
    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 412-419

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 421-430

.. code-block:: Python

    graph = calibrationResult.drawParameterDistributions()
    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.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 431-467

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

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

.. GENERATED FROM PYTHON SOURCE LINES 469-478

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

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

.. GENERATED FROM PYTHON SOURCE LINES 483-484

The standard deviation of the observations.

.. GENERATED FROM PYTHON SOURCE LINES 484-486

.. code-block:: Python

    sigmaH = 0.5  # (m^2)








.. GENERATED FROM PYTHON SOURCE LINES 487-488

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

.. GENERATED FROM PYTHON SOURCE LINES 488-491

.. code-block:: Python

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








.. GENERATED FROM PYTHON SOURCE LINES 492-493

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

.. GENERATED FROM PYTHON SOURCE LINES 493-503

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

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

.. GENERATED FROM PYTHON SOURCE LINES 506-510

.. code-block:: Python

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








.. GENERATED FROM PYTHON SOURCE LINES 511-513

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

.. GENERATED FROM PYTHON SOURCE LINES 513-516

.. code-block:: Python

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








.. GENERATED FROM PYTHON SOURCE LINES 517-521

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 521-524

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

.. code-block:: Python

    graph = calibrationResult.drawObservationsVsInputs()
    graph.setLegendPosition("topleft")
    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 530-532

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

.. GENERATED FROM PYTHON SOURCE LINES 534-537

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

In this case, the fit is satisfactory after calibration.

.. GENERATED FROM PYTHON SOURCE LINES 541-545

.. code-block:: Python

    graph = calibrationResult.drawResiduals()
    graph.setLegendPosition("topleft")
    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 546-548

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

.. GENERATED FROM PYTHON SOURCE LINES 550-552

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

.. GENERATED FROM PYTHON SOURCE LINES 554-557

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

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

.. GENERATED FROM PYTHON SOURCE LINES 559-565

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

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 570-571

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

.. GENERATED FROM PYTHON SOURCE LINES 573-581

.. code-block:: Python

    graph = calibrationResult.drawParameterDistributions()
    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.8)




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

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 591-592

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

.. GENERATED FROM PYTHON SOURCE LINES 592-601

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

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

.. GENERATED FROM PYTHON SOURCE LINES 606-610

.. code-block:: Python

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








.. GENERATED FROM PYTHON SOURCE LINES 611-613

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

.. GENERATED FROM PYTHON SOURCE LINES 613-616

.. code-block:: Python

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








.. GENERATED FROM PYTHON SOURCE LINES 617-621

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

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

.. code-block:: Python

    graph = calibrationResult.drawObservationsVsInputs()
    graph.setLegendPosition("topleft")
    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 630-632

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

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

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

.. GENERATED FROM PYTHON SOURCE LINES 641-649

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

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 654-656

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

.. GENERATED FROM PYTHON SOURCE LINES 656-658

.. code-block:: Python

    distributionPosterior = calibrationResult.getParameterPosterior()








.. GENERATED FROM PYTHON SOURCE LINES 659-660

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

.. GENERATED FROM PYTHON SOURCE LINES 660-666

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

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

.. GENERATED FROM PYTHON SOURCE LINES 670-671

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

.. GENERATED FROM PYTHON SOURCE LINES 671-679

.. code-block:: Python

    graph = calibrationResult.drawParameterDistributions()
    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.8)




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

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 684-685

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

.. GENERATED FROM PYTHON SOURCE LINES 685-694

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

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 711-714

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 714-716

.. code-block:: Python

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





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

 .. code-block:: none

    100




.. GENERATED FROM PYTHON SOURCE LINES 717-719

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

.. GENERATED FROM PYTHON SOURCE LINES 719-733

.. 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, 4.0)},
        legend_kw={"bbox_to_anchor": (1.0, 1.0), "loc": "upper left"},
    )
    plt.subplots_adjust(right=0.8)




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

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

.. GENERATED FROM PYTHON SOURCE LINES 737-738

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

.. GENERATED FROM PYTHON SOURCE LINES 738-749

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

Reset default settings

.. GENERATED FROM PYTHON SOURCE LINES 751-752

.. code-block:: Python

    ot.ResourceMap.Reload()








.. rst-class:: sphx-glr-timing

   **Total running time of the script:** (0 minutes 9.589 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 <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>`