.. only:: html

    .. note::
        :class: sphx-glr-download-link-note

        Click :ref:`here <sphx_glr_download_auto_data_analysis_distribution_fitting_plot_fit_extreme_value_distribution.py>`     to download the full example code
    .. rst-class:: sphx-glr-example-title

    .. _sphx_glr_auto_data_analysis_distribution_fitting_plot_fit_extreme_value_distribution.py:


Fit an extreme value distribution
=================================


.. code-block:: default

    from __future__ import print_function
    import openturns as ot
    import openturns.viewer as viewer
    from matplotlib import pylab as plt

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








Set the random generator seed


.. code-block:: default

    ot.RandomGenerator.SetSeed(0)









The generalized extreme value distribution (GEV)
------------------------------------------------

The :class:`~openturns.GeneralizedExtremeValue` distribution is a family of continuous probability distributions that combine the :class:`~openturns.Gumbel`, :class:`~openturns.Frechet` and :class:`~openturns.WeibullMax` distribution, all extreme value distribution.

In this example we use the associated :class:`~openturns.GeneralizedExtremeValueFactory` to fit sample with extreme values. This factory returns the best model among the Frechet, Gumbel and Weibull estimates according to the BIC criterion.


For experiment purpose we draw a sample from a Gumbel of parameters :math:`\beta = 1.0` and :math:`\gamma = 3.0` and another one from a Frechet with parameters :math:`\beta=1.0, \alpha = 1.0` and :math:`\gamma = 0.0` .
We consider both samples as a unique sample from an unknown extreme distribution to be fitted.


The distributions used :


.. code-block:: default

    myGumbel = ot.Gumbel(1.0, 3.0)
    myFrechet = ot.Frechet(1.0, 1.0, 0.0)









We build our experiment sample of size 2000.


.. code-block:: default

    sample = ot.Sample()
    sampleFrechet = myFrechet.getSample(1000)
    sampleGumbel = myGumbel.getSample(1000)
    sample.add(sampleFrechet)
    sample.add(sampleGumbel)









We fit the sample thanks to the `GeneralizedExtremeValueFactory` :


.. code-block:: default

    myDistribution = ot.GeneralizedExtremeValueFactory().buildAsGeneralizedExtremeValue(
        sample
    )








We can display the parameters of the fitted distribution `myDistribution` :


.. code-block:: default

    print(myDistribution)





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

 Out:

 .. code-block:: none

    GeneralizedExtremeValue(mu=1.79565, sigma=1.54463, xi=0.546359)




We can also get the actual distribution (Weibull, Frechet or Gumbel) with the `getActualDistribution` method :


.. code-block:: default

    print(myDistribution.getActualDistribution())





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

 Out:

 .. code-block:: none

    Frechet(beta = 2.82713, alpha = 1.8303, gamma = -1.03148)




The given sample is then best described by a WeibullMax distribution.

We draw the fitted distribution and an histogram of the data.


.. code-block:: default

    graph = myDistribution.drawPDF()
    graph.add(ot.HistogramFactory().build(sample).drawPDF())
    graph.setColors(["black", "red"])
    graph.setLegends(["GEV fitting", "histogram"])
    graph.setLegendPosition("topright")

    view = viewer.View(graph)
    axes = view.getAxes()
    _ = axes[0].set_xlim(-20.0, 20.0)





.. image:: /auto_data_analysis/distribution_fitting/images/sphx_glr_plot_fit_extreme_value_distribution_001.png
    :alt: plot fit extreme value distribution
    :class: sphx-glr-single-img





For educational purpose we compare different fitting strategies for this sample :

- we use the histogram from the data (in red)
- the GEV fitted distribution (in black)
- the pure Frechet fitted distribution (in green)
- the pure Gumbel fitted distribution (in blue)
- the pure WeibullMax fitted distribution (in dashed orange)



.. code-block:: default

    graph = myDistribution.drawPDF()
    graph.add(ot.HistogramFactory().build(sample).drawPDF())

    distFrechet = ot.FrechetFactory().buildAsFrechet(sample)
    graph.add(distFrechet.drawPDF())

    distGumbel = ot.GumbelFactory().buildAsGumbel(sample)
    graph.add(distGumbel.drawPDF())

    # We change the line style of the WeibullMax because it is expected to be the same 
    # as the GEV one.
    distWeibullMax = ot.WeibullMaxFactory().buildAsWeibullMax(sample)
    curveWeibullMax = distWeibullMax.drawPDF().getDrawable(0)
    curveWeibullMax.setLineStyle("dashed")
    graph.add(curveWeibullMax)

    graph.setColors(["black", "red", "green", "blue", "orange"])
    graph.setLegends(
        [
            "GEV fitting",
            "histogram",
            "Frechet fitting",
            "Gumbel fitting",
            "WeibullMax fitting",
        ]
    )
    graph.setLegendPosition("topright")
    view = viewer.View(graph)
    axes = view.getAxes()  # axes is a matplotlib object
    _ = axes[0].set_xlim(-20.0, 20.0)





.. image:: /auto_data_analysis/distribution_fitting/images/sphx_glr_plot_fit_extreme_value_distribution_002.png
    :alt: plot fit extreme value distribution
    :class: sphx-glr-single-img





As returned by the `getActualDistribution` method the GEV distribution is a WeibullMax.
The :class:`~openturns.GeneralizedExtremeValueFactory` class is a convenient class to fit extreme valued samples without an a priori knowledge of the underlying (at least the closest) extreme distribution.

The Generalized Pareto Distribution (GPD)
-----------------------------------------

In this paragraph we turn to the fitting of a :class:`~openturns.GeneralizedPareto` distribution.
Various estimators are defined in OpenTurns for the GPD factory. Please refer to the :class:`~openturns.GeneralizedParetoFactory` class documentation for more information.
The selection is based on the sample size and compared to the `GeneralizedParetoFactory-SmallSize` key of the :class:`~openturns.ResourceMap` :



.. code-block:: default

    print(ot.ResourceMap.GetAsUnsignedInteger("GeneralizedParetoFactory-SmallSize"))






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

 Out:

 .. code-block:: none

    20




The small sample case
^^^^^^^^^^^^^^^^^^^^^

In this case the default estimator is based on a probability weighted method of moments, with a fallback on the exponential regression method.


.. code-block:: default

    myDist = ot.GeneralizedPareto(1.0, 0.0, 0.0)
    N = 17
    sample = myDist.getSample(N)








We build our experiment sample of size N


.. code-block:: default

    myFittedDist = ot.GeneralizedParetoFactory().buildAsGeneralizedPareto(sample)
    print(myFittedDist)





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

 Out:

 .. code-block:: none

    GeneralizedPareto(sigma = 0.678732, xi=0.0289962, u=0.0498077)




We draw the fitted distribution as well as an histogram to visualize the fitting :


.. code-block:: default

    graph = myFittedDist.drawPDF()
    graph.add(ot.HistogramFactory().build(sample).drawPDF())
    graph.setTitle("Generalized Pareto distribution fitting on a sample")
    graph.setColors(["black", "red"])
    graph.setLegends(["GPD fitting", "histogram"])
    graph.setLegendPosition("topright")

    view = viewer.View(graph)
    axes = view.getAxes()
    _ = axes[0].set_xlim(-1.0, 10.0)





.. image:: /auto_data_analysis/distribution_fitting/images/sphx_glr_plot_fit_extreme_value_distribution_003.png
    :alt: Generalized Pareto distribution fitting on a sample
    :class: sphx-glr-single-img





Large sample
^^^^^^^^^^^^

For a larger sample the estimator is based on an exponential regression method
with a fallback on the probability weighted moments estimator.


.. code-block:: default

    N = 1000
    sample = myDist.getSample(N)








We build our experiment sample of size N


.. code-block:: default

    myFittedDist = ot.GeneralizedParetoFactory().buildAsGeneralizedPareto(sample)
    print(myFittedDist)





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

 Out:

 .. code-block:: none

    GeneralizedPareto(sigma = 0.971553, xi=-0.000639593, u=0.000103683)




We draw the fitted distribution as well as an histogram to visualize the fitting :


.. code-block:: default

    graph = myFittedDist.drawPDF()
    graph.add(ot.HistogramFactory().build(sample).drawPDF())
    graph.setTitle("Generalized Pareto distribution fitting on a sample")
    graph.setColors(["black", "red"])
    graph.setLegends(["GPD fitting", "histogram"])
    graph.setLegendPosition("topright")

    view = viewer.View(graph)
    axes = view.getAxes()
    _ = axes[0].set_xlim(-1.0, 10.0)


    plt.show()



.. image:: /auto_data_analysis/distribution_fitting/images/sphx_glr_plot_fit_extreme_value_distribution_004.png
    :alt: Generalized Pareto distribution fitting on a sample
    :class: sphx-glr-single-img






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

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


.. _sphx_glr_download_auto_data_analysis_distribution_fitting_plot_fit_extreme_value_distribution.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_fit_extreme_value_distribution.py <plot_fit_extreme_value_distribution.py>`



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

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


.. only:: html

 .. rst-class:: sphx-glr-signature

    `Gallery generated by Sphinx-Gallery <https://sphinx-gallery.github.io>`_