.. only:: html .. note:: :class: sphx-glr-download-link-note Click :ref:`here ` 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 ================================= We present a calibration study of the logistic growth model defined :ref:`here `. Analysis of the data -------------------- .. code-block:: default 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) We load the logistic model from the usecases module : .. code-block:: default from openturns.usecases import logistic_model as logistic_model lm = logistic_model.LogisticModel() The data is based on 22 dates from 1790 to 2000. The observation points are stored in the `data` field : .. code-block:: default observedSample = lm.data .. code-block:: default nbobs = observedSample.getSize() nbobs .. rst-class:: sphx-glr-script-out Out: .. code-block:: none 22 .. code-block:: default timeObservations = observedSample[:,0] timeObservations[0:5] .. raw:: html
v0
01790
11800
21810
31820
41830


.. code-block:: default populationObservations = observedSample[:,1] populationObservations[0:5] .. raw:: html
v1
03.9
15.3
27.2
39.6
413


.. 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:: /auto_calibration/least_squares_and_gaussian_calibration/images/sphx_glr_plot_calibration_logistic_001.png :alt: plot calibration logistic :class: sphx-glr-single-img 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`. .. code-block:: default nbdates = 22 def logisticModel(X): t = [X[i] for i in range(nbdates)] a = X[22] c = X[23] t0 = 1790. y0 = 3.9e6 b = np.exp(c) y = ot.Point(nbdates) for i in range(nbdates): y[i] = a*y0/(b*y0+(a-b*y0)*np.exp(-a*(t[i]-t0))) z = y/1.e6 # Convert into millions return z .. code-block:: default logisticModelPy = ot.PythonFunction(24, nbdates, logisticModel) The reference values of the parameters. Because :math:`b` is so small with respect to :math:`a`, we use the logarithm. .. code-block:: default np.log(1.5587e-10) .. rst-class:: sphx-glr-script-out Out: .. code-block:: none -22.581998789427587 .. code-block:: default a=0.03134 c=-22.58 thetaPrior = [a,c] .. code-block:: default logisticParametric = ot.ParametricFunction(logisticModelPy,[22,23],thetaPrior) Check that we can evaluate the parametric function. .. code-block:: default populationPredicted = logisticParametric(timeObservations.asPoint()) populationPredicted .. raw:: html

[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



.. 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:: /auto_calibration/least_squares_and_gaussian_calibration/images/sphx_glr_plot_calibration_logistic_002.png :alt: plot calibration logistic :class: sphx-glr-single-img We see that the fit is not good: the observations continue to grow after 1950, while the growth of the prediction seem to fade. Calibration with linear least squares ------------------------------------- .. code-block:: default timeObservationsVector = ot.Sample([[timeObservations[i,0] for i in range(nbobs)]]) timeObservationsVector[0:10] .. raw:: html
v0v1v2v3v4v5v6v7v8v9v10v11v12v13v14v15v16v17v18v19v20v21
01790180018101820183018401850186018701880189019001910192019301940195019601970198019902000


.. code-block:: default populationObservationsVector = ot.Sample([[populationObservations[i, 0] for i in range(nbobs)]]) populationObservationsVector[0:10] .. raw:: html
v0v1v2v3v4v5v6v7v8v9v10v11v12v13v14v15v16v17v18v19v20v21
03.95.37.29.6131723313950627692106123132151179203221250281


The reference values of the parameters. .. code-block:: default a=0.03134 c=-22.58 thetaPrior = ot.Point([a,c]) .. code-block:: default logisticParametric = ot.ParametricFunction(logisticModelPy,[22,23],thetaPrior) Check that we can evaluate the parametric function. .. code-block:: default populationPredicted = logisticParametric(timeObservationsVector) populationPredicted[0:10] .. raw:: html
y0y1y2y3y4y5y6y7y8y9y10y11y12y13y14y15y16y17y18y19y20y21
03.95.2975717.1776949.69197713.0276917.4068223.0769130.288739.2560650.0976762.7690777.006392.31103108.0009123.3223137.5899150.3003161.1843170.193177.4422183.1443187.5496


Calibration ------------ .. code-block:: default algo = ot.LinearLeastSquaresCalibration(logisticParametric, timeObservationsVector, populationObservationsVector, thetaPrior) .. code-block:: default algo.run() .. code-block:: default calibrationResult = algo.getResult() .. code-block:: default thetaMAP = calibrationResult.getParameterMAP() thetaMAP .. raw:: html

[0.0265958,-23.1714]



.. code-block:: default thetaPosterior = calibrationResult.getParameterPosterior() thetaPosterior.computeBilateralConfidenceIntervalWithMarginalProbability(0.95)[0] .. raw:: html

[0.0246465, 0.028545]
[-23.3182, -23.0247]



transpose samples to interpret several observations instead of several input/outputs as it is a field model .. 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) .. code-block:: default graph = calibrationResult.drawObservationsVsInputs() view = viewer.View(graph) .. image:: /auto_calibration/least_squares_and_gaussian_calibration/images/sphx_glr_plot_calibration_logistic_003.png :alt: plot calibration logistic :class: sphx-glr-single-img .. code-block:: default graph = calibrationResult.drawObservationsVsInputs() view = viewer.View(graph) .. image:: /auto_calibration/least_squares_and_gaussian_calibration/images/sphx_glr_plot_calibration_logistic_004.png :alt: plot calibration logistic :class: sphx-glr-single-img .. code-block:: default graph = calibrationResult.drawObservationsVsPredictions() view = viewer.View(graph) .. image:: /auto_calibration/least_squares_and_gaussian_calibration/images/sphx_glr_plot_calibration_logistic_005.png :alt: plot calibration logistic :class: sphx-glr-single-img .. code-block:: default graph = calibrationResult.drawResiduals() view = viewer.View(graph) .. image:: /auto_calibration/least_squares_and_gaussian_calibration/images/sphx_glr_plot_calibration_logistic_006.png :alt: plot calibration logistic :class: sphx-glr-single-img .. code-block:: default graph = calibrationResult.drawParameterDistributions() view = viewer.View(graph) plt.show() .. image:: /auto_calibration/least_squares_and_gaussian_calibration/images/sphx_glr_plot_calibration_logistic_007.png :alt: plot calibration logistic :class: sphx-glr-single-img .. rst-class:: sphx-glr-timing **Total running time of the script:** ( 0 minutes 1.315 seconds) .. _sphx_glr_download_auto_calibration_least_squares_and_gaussian_calibration_plot_calibration_logistic.py: .. only :: html .. container:: sphx-glr-footer :class: sphx-glr-footer-example .. container:: sphx-glr-download sphx-glr-download-python :download:`Download Python source code: plot_calibration_logistic.py ` .. container:: sphx-glr-download sphx-glr-download-jupyter :download:`Download Jupyter notebook: plot_calibration_logistic.ipynb ` .. only:: html .. rst-class:: sphx-glr-signature `Gallery generated by Sphinx-Gallery `_