.. 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_logistic.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_logistic.py>`
        to download the full example code

.. rst-class:: sphx-glr-example-title

.. _sphx_glr_auto_calibration_least_squares_and_gaussian_calibration_plot_calibration_logistic.py:


Calibration of the logistic model
=================================

.. GENERATED FROM PYTHON SOURCE LINES 6-8

We present a calibration study of the logistic growth model defined :ref:`here <use-case-logistic>`.
In this example, we calibrate the parameters of a model which predicts the dynamics of the size of a population. This shows how you can calibrate a model which predicts a time dependent output. You need to view the output time series as a vector.

.. GENERATED FROM PYTHON SOURCE LINES 10-12

Analysis of the data
--------------------

.. GENERATED FROM PYTHON SOURCE LINES 14-23

.. code-block:: default

    from openturns.usecases import logistic_model
    import openturns as ot
    import numpy as np
    import openturns.viewer as viewer
    from matplotlib import pylab as plt

    ot.Log.Show(ot.Log.NONE)









.. GENERATED FROM PYTHON SOURCE LINES 24-25

We load the logistic model from the usecases module :

.. GENERATED FROM PYTHON SOURCE LINES 25-27

.. code-block:: default

    lm = logistic_model.LogisticModel()








.. GENERATED FROM PYTHON SOURCE LINES 28-29

The data is based on 22 dates from 1790 to 2000. The observation points are stored in the `data` field :

.. GENERATED FROM PYTHON SOURCE LINES 29-31

.. code-block:: default

    observedSample = lm.data








.. GENERATED FROM PYTHON SOURCE LINES 32-35

.. code-block:: default

    nbobs = observedSample.getSize()
    nbobs





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

 .. code-block:: none


    22



.. GENERATED FROM PYTHON SOURCE LINES 36-39

.. code-block:: default

    timeObservations = observedSample[:, 0]
    timeObservations[0:5]






.. raw:: html

    <div class="output_subarea output_html rendered_html output_result">
    <TABLE><TR><TD></TD><TH>v0</TH></TR>
    <TR><TH>0</TH><TD>1790</TD></TR>
    <TR><TH>1</TH><TD>1800</TD></TR>
    <TR><TH>2</TH><TD>1810</TD></TR>
    <TR><TH>3</TH><TD>1820</TD></TR>
    <TR><TH>4</TH><TD>1830</TD></TR>
    </TABLE>
    </div>
    <br />
    <br />

.. GENERATED FROM PYTHON SOURCE LINES 40-43

.. code-block:: default

    populationObservations = observedSample[:, 1]
    populationObservations[0:5]






.. raw:: html

    <div class="output_subarea output_html rendered_html output_result">
    <TABLE><TR><TD></TD><TH>v1</TH></TR>
    <TR><TH>0</TH><TD>3.9</TD></TR>
    <TR><TH>1</TH><TD>5.3</TD></TR>
    <TR><TH>2</TH><TD>7.2</TD></TR>
    <TR><TH>3</TH><TD>9.6</TD></TR>
    <TR><TH>4</TH><TD>13</TD></TR>
    </TABLE>
    </div>
    <br />
    <br />

.. GENERATED FROM PYTHON SOURCE LINES 44-50

.. code-block:: default

    graph = ot.Graph("", "Time (years)", "Population (Millions)", True, "topleft")
    cloud = ot.Cloud(timeObservations, populationObservations)
    cloud.setLegend("Observations")
    graph.add(cloud)
    view = viewer.View(graph)




.. image-sg:: /auto_calibration/least_squares_and_gaussian_calibration/images/sphx_glr_plot_calibration_logistic_001.png
   :alt: plot calibration logistic
   :srcset: /auto_calibration/least_squares_and_gaussian_calibration/images/sphx_glr_plot_calibration_logistic_001.png
   :class: sphx-glr-single-img





.. GENERATED FROM PYTHON SOURCE LINES 51-52

We consider the times and populations as dimension 22 vectors. The `logisticModel` function takes a dimension 24 vector as input and returns a dimension 22 vector. The first 22 components of the input vector contains the dates and the remaining 2 components are :math:`a` and :math:`b`.

.. GENERATED FROM PYTHON SOURCE LINES 54-71

.. code-block:: default

    nbdates = 22


    def logisticModel(X):
        t = [X[i] for i in range(nbdates)]
        a = X[22]
        c = X[23]
        t0 = 1790.0
        y0 = 3.9e6
        b = np.exp(c)
        y = [0.0] * nbdates
        for i in range(nbdates):
            y[i] = a * y0 / (b * y0 + (a - b * y0) * np.exp(-a * (t[i] - t0)))
        z = [yi / 1.0e6 for yi in y]  # Convert into millions
        return z









.. GENERATED FROM PYTHON SOURCE LINES 72-74

.. code-block:: default

    logisticModelPy = ot.PythonFunction(24, nbdates, logisticModel)








.. GENERATED FROM PYTHON SOURCE LINES 75-78

The reference values of the parameters.

Because :math:`b` is so small with respect to :math:`a`, we use the logarithm.

.. GENERATED FROM PYTHON SOURCE LINES 80-82

.. code-block:: default

    np.log(1.5587e-10)





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

 .. code-block:: none


    -22.581998789427587



.. GENERATED FROM PYTHON SOURCE LINES 83-87

.. code-block:: default

    a = 0.03134
    c = -22.58
    thetaPrior = [a, c]








.. GENERATED FROM PYTHON SOURCE LINES 88-90

.. code-block:: default

    logisticParametric = ot.ParametricFunction(logisticModelPy, [22, 23], thetaPrior)








.. GENERATED FROM PYTHON SOURCE LINES 91-92

Check that we can evaluate the parametric function.

.. GENERATED FROM PYTHON SOURCE LINES 94-97

.. code-block:: default

    populationPredicted = logisticParametric(timeObservations.asPoint())
    populationPredicted






.. raw:: html

    <div class="output_subarea output_html rendered_html output_result">
    <p>[3.9,5.29757,7.17769,9.69198,13.0277,17.4068,23.0769,30.2887,39.2561,50.0977,62.7691,77.0063,92.311,108.001,123.322,137.59,150.3,161.184,170.193,177.442,183.144,187.55]#22</p>
    </div>
    <br />
    <br />

.. GENERATED FROM PYTHON SOURCE LINES 98-111

.. code-block:: default

    graph = ot.Graph("", "Time (years)", "Population (Millions)", True, "topleft")
    # Observations
    cloud = ot.Cloud(timeObservations, populationObservations)
    cloud.setLegend("Observations")
    cloud.setColor("blue")
    graph.add(cloud)
    # Predictions
    cloud = ot.Cloud(timeObservations.asPoint(), populationPredicted)
    cloud.setLegend("Predictions")
    cloud.setColor("green")
    graph.add(cloud)
    view = viewer.View(graph)




.. image-sg:: /auto_calibration/least_squares_and_gaussian_calibration/images/sphx_glr_plot_calibration_logistic_002.png
   :alt: plot calibration logistic
   :srcset: /auto_calibration/least_squares_and_gaussian_calibration/images/sphx_glr_plot_calibration_logistic_002.png
   :class: sphx-glr-single-img





.. GENERATED FROM PYTHON SOURCE LINES 112-113

We see that the fit is not good: the observations continue to grow after 1950, while the growth of the prediction seem to fade.

.. GENERATED FROM PYTHON SOURCE LINES 115-117

Calibration with linear least squares
-------------------------------------

.. GENERATED FROM PYTHON SOURCE LINES 119-122

.. code-block:: default

    timeObservationsVector = ot.Sample([[timeObservations[i, 0] for i in range(nbobs)]])
    timeObservationsVector[0:10]






.. raw:: html

    <div class="output_subarea output_html rendered_html output_result">
    <TABLE><TR><TD></TD><TH>v0</TH><TH>v1</TH><TH>v2</TH><TH>v3</TH><TH>v4</TH><TH>v5</TH><TH>v6</TH><TH>v7</TH><TH>v8</TH><TH>v9</TH><TH>v10</TH><TH>v11</TH><TH>v12</TH><TH>v13</TH><TH>v14</TH><TH>v15</TH><TH>v16</TH><TH>v17</TH><TH>v18</TH><TH>v19</TH><TH>v20</TH><TH>v21</TH></TR>
    <TR><TH>0</TH><TD>1790</TD><TD>1800</TD><TD>1810</TD><TD>1820</TD><TD>1830</TD><TD>1840</TD><TD>1850</TD><TD>1860</TD><TD>1870</TD><TD>1880</TD><TD>1890</TD><TD>1900</TD><TD>1910</TD><TD>1920</TD><TD>1930</TD><TD>1940</TD><TD>1950</TD><TD>1960</TD><TD>1970</TD><TD>1980</TD><TD>1990</TD><TD>2000</TD></TR>
    </TABLE>
    </div>
    <br />
    <br />

.. GENERATED FROM PYTHON SOURCE LINES 123-128

.. code-block:: default

    populationObservationsVector = ot.Sample(
        [[populationObservations[i, 0] for i in range(nbobs)]]
    )
    populationObservationsVector[0:10]






.. raw:: html

    <div class="output_subarea output_html rendered_html output_result">
    <TABLE><TR><TD></TD><TH>v0</TH><TH>v1</TH><TH>v2</TH><TH>v3</TH><TH>v4</TH><TH>v5</TH><TH>v6</TH><TH>v7</TH><TH>v8</TH><TH>v9</TH><TH>v10</TH><TH>v11</TH><TH>v12</TH><TH>v13</TH><TH>v14</TH><TH>v15</TH><TH>v16</TH><TH>v17</TH><TH>v18</TH><TH>v19</TH><TH>v20</TH><TH>v21</TH></TR>
    <TR><TH>0</TH><TD>3.9</TD><TD>5.3</TD><TD>7.2</TD><TD>9.6</TD><TD>13</TD><TD>17</TD><TD>23</TD><TD>31</TD><TD>39</TD><TD>50</TD><TD>62</TD><TD>76</TD><TD>92</TD><TD>106</TD><TD>123</TD><TD>132</TD><TD>151</TD><TD>179</TD><TD>203</TD><TD>221</TD><TD>250</TD><TD>281</TD></TR>
    </TABLE>
    </div>
    <br />
    <br />

.. GENERATED FROM PYTHON SOURCE LINES 129-130

The reference values of the parameters.

.. GENERATED FROM PYTHON SOURCE LINES 132-137

.. code-block:: default

    a = 0.03134
    c = -22.58
    thetaPrior = [a, c]









.. GENERATED FROM PYTHON SOURCE LINES 138-140

.. code-block:: default

    logisticParametric = ot.ParametricFunction(logisticModelPy, [22, 23], thetaPrior)








.. GENERATED FROM PYTHON SOURCE LINES 141-142

Check that we can evaluate the parametric function.

.. GENERATED FROM PYTHON SOURCE LINES 144-147

.. code-block:: default

    populationPredicted = logisticParametric(timeObservationsVector)
    populationPredicted[0:10]






.. raw:: html

    <div class="output_subarea output_html rendered_html output_result">
    <TABLE><TR><TD></TD><TH>y0</TH><TH>y1</TH><TH>y2</TH><TH>y3</TH><TH>y4</TH><TH>y5</TH><TH>y6</TH><TH>y7</TH><TH>y8</TH><TH>y9</TH><TH>y10</TH><TH>y11</TH><TH>y12</TH><TH>y13</TH><TH>y14</TH><TH>y15</TH><TH>y16</TH><TH>y17</TH><TH>y18</TH><TH>y19</TH><TH>y20</TH><TH>y21</TH></TR>
    <TR><TH>0</TH><TD>3.9</TD><TD>5.297571</TD><TD>7.177694</TD><TD>9.691977</TD><TD>13.02769</TD><TD>17.40682</TD><TD>23.07691</TD><TD>30.2887</TD><TD>39.25606</TD><TD>50.09767</TD><TD>62.76907</TD><TD>77.0063</TD><TD>92.31103</TD><TD>108.0009</TD><TD>123.3223</TD><TD>137.5899</TD><TD>150.3003</TD><TD>161.1843</TD><TD>170.193</TD><TD>177.4422</TD><TD>183.1443</TD><TD>187.5496</TD></TR>
    </TABLE>
    </div>
    <br />
    <br />

.. GENERATED FROM PYTHON SOURCE LINES 148-150

Calibration
------------

.. GENERATED FROM PYTHON SOURCE LINES 152-156

.. code-block:: default

    algo = ot.LinearLeastSquaresCalibration(
        logisticParametric, timeObservationsVector, populationObservationsVector, thetaPrior
    )








.. GENERATED FROM PYTHON SOURCE LINES 157-159

.. code-block:: default

    algo.run()








.. GENERATED FROM PYTHON SOURCE LINES 160-162

.. code-block:: default

    calibrationResult = algo.getResult()








.. GENERATED FROM PYTHON SOURCE LINES 163-166

.. code-block:: default

    thetaMAP = calibrationResult.getParameterMAP()
    thetaMAP






.. raw:: html

    <div class="output_subarea output_html rendered_html output_result">
    <p>[0.0265958,-23.1714]</p>
    </div>
    <br />
    <br />

.. GENERATED FROM PYTHON SOURCE LINES 167-170

.. code-block:: default

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






.. raw:: html

    <div class="output_subarea output_html rendered_html output_result">
    <p>[0.0246465, 0.028545]<br>
    [-23.3182, -23.0247]</p>
    </div>
    <br />
    <br />

.. GENERATED FROM PYTHON SOURCE LINES 171-172

transpose samples to interpret several observations instead of several input/outputs as it is a field model

.. GENERATED FROM PYTHON SOURCE LINES 172-187

.. code-block:: default

    if calibrationResult.getInputObservations().getSize() == 1:
        calibrationResult.setInputObservations(
            [timeObservations[i] for i in range(nbdates)]
        )
        calibrationResult.setOutputObservations(
            [populationObservations[i] for i in range(nbdates)]
        )
        outputAtPrior = [
            [calibrationResult.getOutputAtPriorMean()[0, i]] for i in range(nbdates)
        ]
        outputAtPosterior = [
            [calibrationResult.getOutputAtPosteriorMean()[0, i]] for i in range(nbdates)
        ]
        calibrationResult.setOutputAtPriorAndPosteriorMean(outputAtPrior, outputAtPosterior)








.. GENERATED FROM PYTHON SOURCE LINES 188-191

.. code-block:: default

    graph = calibrationResult.drawObservationsVsInputs()
    view = viewer.View(graph)




.. image-sg:: /auto_calibration/least_squares_and_gaussian_calibration/images/sphx_glr_plot_calibration_logistic_003.png
   :alt: plot calibration logistic
   :srcset: /auto_calibration/least_squares_and_gaussian_calibration/images/sphx_glr_plot_calibration_logistic_003.png
   :class: sphx-glr-single-img





.. GENERATED FROM PYTHON SOURCE LINES 192-195

.. code-block:: default

    graph = calibrationResult.drawObservationsVsInputs()
    view = viewer.View(graph)




.. image-sg:: /auto_calibration/least_squares_and_gaussian_calibration/images/sphx_glr_plot_calibration_logistic_004.png
   :alt: plot calibration logistic
   :srcset: /auto_calibration/least_squares_and_gaussian_calibration/images/sphx_glr_plot_calibration_logistic_004.png
   :class: sphx-glr-single-img





.. GENERATED FROM PYTHON SOURCE LINES 196-199

.. code-block:: default

    graph = calibrationResult.drawObservationsVsPredictions()
    view = viewer.View(graph)




.. image-sg:: /auto_calibration/least_squares_and_gaussian_calibration/images/sphx_glr_plot_calibration_logistic_005.png
   :alt: plot calibration logistic
   :srcset: /auto_calibration/least_squares_and_gaussian_calibration/images/sphx_glr_plot_calibration_logistic_005.png
   :class: sphx-glr-single-img





.. GENERATED FROM PYTHON SOURCE LINES 200-203

.. code-block:: default

    graph = calibrationResult.drawResiduals()
    view = viewer.View(graph)




.. image-sg:: /auto_calibration/least_squares_and_gaussian_calibration/images/sphx_glr_plot_calibration_logistic_006.png
   :alt: , Residual analysis
   :srcset: /auto_calibration/least_squares_and_gaussian_calibration/images/sphx_glr_plot_calibration_logistic_006.png
   :class: sphx-glr-single-img





.. GENERATED FROM PYTHON SOURCE LINES 204-208

.. code-block:: default

    graph = calibrationResult.drawParameterDistributions()
    view = viewer.View(graph)

    plt.show()



.. image-sg:: /auto_calibration/least_squares_and_gaussian_calibration/images/sphx_glr_plot_calibration_logistic_007.png
   :alt: plot calibration logistic
   :srcset: /auto_calibration/least_squares_and_gaussian_calibration/images/sphx_glr_plot_calibration_logistic_007.png
   :class: sphx-glr-single-img






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

   **Total running time of the script:** ( 0 minutes  0.851 seconds)


.. _sphx_glr_download_auto_calibration_least_squares_and_gaussian_calibration_plot_calibration_logistic.py:

.. only:: html

  .. container:: sphx-glr-footer sphx-glr-footer-example




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

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

    .. container:: sphx-glr-download sphx-glr-download-jupyter

      :download:`Download Jupyter notebook: plot_calibration_logistic.ipynb <plot_calibration_logistic.ipynb>`