`.
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 |
0 | 1790 |
1 | 1800 |
2 | 1810 |
3 | 1820 |
4 | 1830 |
.. code-block:: default
populationObservations = observedSample[:,1]
populationObservations[0:5]
.. raw:: html
| v1 |
0 | 3.9 |
1 | 5.3 |
2 | 7.2 |
3 | 9.6 |
4 | 13 |
.. 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
| v0 | v1 | v2 | v3 | v4 | v5 | v6 | v7 | v8 | v9 | v10 | v11 | v12 | v13 | v14 | v15 | v16 | v17 | v18 | v19 | v20 | v21 |
0 | 1790 | 1800 | 1810 | 1820 | 1830 | 1840 | 1850 | 1860 | 1870 | 1880 | 1890 | 1900 | 1910 | 1920 | 1930 | 1940 | 1950 | 1960 | 1970 | 1980 | 1990 | 2000 |
.. code-block:: default
populationObservationsVector = ot.Sample([[populationObservations[i, 0] for i in range(nbobs)]])
populationObservationsVector[0:10]
.. raw:: html
| v0 | v1 | v2 | v3 | v4 | v5 | v6 | v7 | v8 | v9 | v10 | v11 | v12 | v13 | v14 | v15 | v16 | v17 | v18 | v19 | v20 | v21 |
0 | 3.9 | 5.3 | 7.2 | 9.6 | 13 | 17 | 23 | 31 | 39 | 50 | 62 | 76 | 92 | 106 | 123 | 132 | 151 | 179 | 203 | 221 | 250 | 281 |
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
| y0 | y1 | y2 | y3 | y4 | y5 | y6 | y7 | y8 | y9 | y10 | y11 | y12 | y13 | y14 | y15 | y16 | y17 | y18 | y19 | y20 | y21 |
0 | 3.9 | 5.297571 | 7.177694 | 9.691977 | 13.02769 | 17.40682 | 23.07691 | 30.2887 | 39.25606 | 50.09767 | 62.76907 | 77.0063 | 92.31103 | 108.0009 | 123.3223 | 137.5899 | 150.3003 | 161.1843 | 170.193 | 177.4422 | 183.1443 | 187.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 `_