.. DO NOT EDIT.
.. THIS FILE WAS AUTOMATICALLY GENERATED BY SPHINX-GALLERY.
.. TO MAKE CHANGES, EDIT THE SOURCE PYTHON FILE:
.. "auto_probabilistic_modeling/distributions/plot_minimum_volume_level_sets.py"
.. LINE NUMBERS ARE GIVEN BELOW.

.. only:: html

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

        :ref:`Go to the end <sphx_glr_download_auto_probabilistic_modeling_distributions_plot_minimum_volume_level_sets.py>`
        to download the full example code

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

.. _sphx_glr_auto_probabilistic_modeling_distributions_plot_minimum_volume_level_sets.py:


Draw minimum volume level sets
==============================

.. GENERATED FROM PYTHON SOURCE LINES 6-12

.. code-block:: Python

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

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








.. GENERATED FROM PYTHON SOURCE LINES 13-17

Draw minimum volume level set in 1D
-----------------------------------

In this paragraph, we compute the minimum volume level set of a univariate distribution.

.. GENERATED FROM PYTHON SOURCE LINES 20-22

With a Normal, minimum volume LevelSet
^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^

.. GENERATED FROM PYTHON SOURCE LINES 24-26

.. code-block:: Python

    n = ot.Normal()








.. GENERATED FROM PYTHON SOURCE LINES 27-30

.. code-block:: Python

    graph = n.drawPDF()
    view = viewer.View(graph)




.. image-sg:: /auto_probabilistic_modeling/distributions/images/sphx_glr_plot_minimum_volume_level_sets_001.png
   :alt: plot minimum volume level sets
   :srcset: /auto_probabilistic_modeling/distributions/images/sphx_glr_plot_minimum_volume_level_sets_001.png
   :class: sphx-glr-single-img





.. GENERATED FROM PYTHON SOURCE LINES 31-34

We want to compute the minimum volume LevelSet which contains `alpha`=90% of the distribution.
The `threshold` is the value of the PDF corresponding the `alpha`-probability:
the points contained in the LevelSet have a PDF value lower or equal to this threshold.

.. GENERATED FROM PYTHON SOURCE LINES 36-40

.. code-block:: Python

    alpha = 0.9
    levelSet, threshold = n.computeMinimumVolumeLevelSetWithThreshold(alpha)
    threshold





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

 .. code-block:: none


    0.10313564037537133



.. GENERATED FROM PYTHON SOURCE LINES 41-42

The `LevelSet` has a `contains` method. Obviously, the point 0 is in the LevelSet.

.. GENERATED FROM PYTHON SOURCE LINES 44-47

.. code-block:: Python

    levelSet.contains([0.0])






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

 .. code-block:: none


    True



.. GENERATED FROM PYTHON SOURCE LINES 48-68

.. code-block:: Python

    def computeSampleInLevelSet(distribution, levelSet, sampleSize=1000):
        """
        Generate a sample from given distribution.
        Extract the sub-sample which is contained in the levelSet.
        """
        sample = distribution.getSample(sampleSize)
        dim = distribution.getDimension()
        # Get the list of points in the LevelSet.
        inLevelSet = []
        for x in sample:
            if levelSet.contains(x):
                inLevelSet.append(x)
        # Extract the sub-sample of the points in the LevelSet
        numberOfPointsInLevelSet = len(inLevelSet)
        inLevelSetSample = ot.Sample(numberOfPointsInLevelSet, dim)
        for i in range(numberOfPointsInLevelSet):
            inLevelSetSample[i] = inLevelSet[i]
        return inLevelSetSample









.. GENERATED FROM PYTHON SOURCE LINES 69-80

.. code-block:: Python

    def from1Dto2Dsample(oldSample):
        """
        Create a 2D sample from a 1D sample with zero ordinate (for the graph).
        """
        size = oldSample.getSize()
        newSample = ot.Sample(size, 2)
        for i in range(size):
            newSample[i, 0] = oldSample[i, 0]
        return newSample









.. GENERATED FROM PYTHON SOURCE LINES 81-97

.. code-block:: Python

    def drawLevelSet1D(distribution, levelSet, alpha, threshold, sampleSize=100):
        """
        Draw a 1D sample included in a given levelSet.
        The sample is generated from the distribution.
        """
        inLevelSample = computeSampleInLevelSet(distribution, levelSet, sampleSize)
        cloudSample = from1Dto2Dsample(inLevelSample)
        graph = distribution.drawPDF()
        mycloud = ot.Cloud(cloudSample)
        graph.add(mycloud)
        graph.setTitle(
            "%.2f%% of the distribution, sample size = %d, " % (100 * alpha, sampleSize)
        )
        return graph









.. GENERATED FROM PYTHON SOURCE LINES 98-101

.. code-block:: Python

    graph = drawLevelSet1D(n, levelSet, alpha, threshold)
    view = viewer.View(graph)




.. image-sg:: /auto_probabilistic_modeling/distributions/images/sphx_glr_plot_minimum_volume_level_sets_002.png
   :alt: 90.00% of the distribution, sample size = 100, 
   :srcset: /auto_probabilistic_modeling/distributions/images/sphx_glr_plot_minimum_volume_level_sets_002.png
   :class: sphx-glr-single-img





.. GENERATED FROM PYTHON SOURCE LINES 102-104

With a Normal, minimum volume Interval
^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^

.. GENERATED FROM PYTHON SOURCE LINES 106-110

.. code-block:: Python

    interval = n.computeMinimumVolumeInterval(alpha)
    interval







.. raw:: html

    <div class="output_subarea output_html rendered_html output_result">
    class=Interval name=Unnamed dimension=1 lower bound=class=Point name=Unnamed dimension=1 values=[-1.64485] upper bound=class=Point name=Unnamed dimension=1 values=[1.64485] finite lower bound=[1] finite upper bound=[1]
    </div>
    <br />
    <br />

.. GENERATED FROM PYTHON SOURCE LINES 111-129

.. code-block:: Python

    def drawPDFAndInterval1D(distribution, interval, alpha):
        """
        Draw the PDF of the distribution and the lower and upper bounds of an interval.
        """
        xmin = interval.getLowerBound()[0]
        xmax = interval.getUpperBound()[0]
        graph = distribution.drawPDF()
        yvalue = distribution.computePDF(xmin)
        curve = ot.Curve([[xmin, 0.0], [xmin, yvalue], [xmax, yvalue], [xmax, 0.0]])
        curve.setColor("black")
        graph.add(curve)
        graph.setTitle(
            "%.2f%% of the distribution, lower bound = %.3f, upper bound = %.3f"
            % (100 * alpha, xmin, xmax)
        )
        return graph









.. GENERATED FROM PYTHON SOURCE LINES 130-131

The `computeMinimumVolumeInterval` returns an `Interval`.

.. GENERATED FROM PYTHON SOURCE LINES 133-136

.. code-block:: Python

    graph = drawPDFAndInterval1D(n, interval, alpha)
    view = viewer.View(graph)




.. image-sg:: /auto_probabilistic_modeling/distributions/images/sphx_glr_plot_minimum_volume_level_sets_003.png
   :alt: 90.00% of the distribution, lower bound = -1.645, upper bound = 1.645
   :srcset: /auto_probabilistic_modeling/distributions/images/sphx_glr_plot_minimum_volume_level_sets_003.png
   :class: sphx-glr-single-img





.. GENERATED FROM PYTHON SOURCE LINES 137-139

With a Mixture, minimum volume LevelSet
^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^

.. GENERATED FROM PYTHON SOURCE LINES 141-143

.. code-block:: Python

    m = ot.Mixture([ot.Normal(-5.0, 1.0), ot.Normal(5.0, 1.0)], [0.2, 0.8])








.. GENERATED FROM PYTHON SOURCE LINES 144-147

.. code-block:: Python

    graph = m.drawPDF()
    view = viewer.View(graph)




.. image-sg:: /auto_probabilistic_modeling/distributions/images/sphx_glr_plot_minimum_volume_level_sets_004.png
   :alt: plot minimum volume level sets
   :srcset: /auto_probabilistic_modeling/distributions/images/sphx_glr_plot_minimum_volume_level_sets_004.png
   :class: sphx-glr-single-img





.. GENERATED FROM PYTHON SOURCE LINES 148-152

.. code-block:: Python

    alpha = 0.9
    levelSet, threshold = m.computeMinimumVolumeLevelSetWithThreshold(alpha)
    threshold





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

 .. code-block:: none


    0.04667473141153258



.. GENERATED FROM PYTHON SOURCE LINES 153-154

The interesting point is that a `LevelSet` may be non-contiguous. In the current mixture example, this is not an interval.

.. GENERATED FROM PYTHON SOURCE LINES 156-159

.. code-block:: Python

    graph = drawLevelSet1D(m, levelSet, alpha, threshold, 1000)
    view = viewer.View(graph)




.. image-sg:: /auto_probabilistic_modeling/distributions/images/sphx_glr_plot_minimum_volume_level_sets_005.png
   :alt: 90.00% of the distribution, sample size = 1000, 
   :srcset: /auto_probabilistic_modeling/distributions/images/sphx_glr_plot_minimum_volume_level_sets_005.png
   :class: sphx-glr-single-img





.. GENERATED FROM PYTHON SOURCE LINES 160-162

With a Mixture, minimum volume Interval
^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^

.. GENERATED FROM PYTHON SOURCE LINES 164-167

.. code-block:: Python

    interval = m.computeMinimumVolumeInterval(alpha)
    interval






.. raw:: html

    <div class="output_subarea output_html rendered_html output_result">
    class=Interval name=Unnamed dimension=1 lower bound=class=Point name=Unnamed dimension=1 values=[-5.44003] upper bound=class=Point name=Unnamed dimension=1 values=[6.72227] finite lower bound=[1] finite upper bound=[1]
    </div>
    <br />
    <br />

.. GENERATED FROM PYTHON SOURCE LINES 168-169

The `computeMinimumVolumeInterval` returns an `Interval`. The bounds of this interval are different from the previous `LevelSet`.

.. GENERATED FROM PYTHON SOURCE LINES 171-175

.. code-block:: Python

    graph = drawPDFAndInterval1D(m, interval, alpha)
    view = viewer.View(graph)





.. image-sg:: /auto_probabilistic_modeling/distributions/images/sphx_glr_plot_minimum_volume_level_sets_006.png
   :alt: 90.00% of the distribution, lower bound = -5.440, upper bound = 6.722
   :srcset: /auto_probabilistic_modeling/distributions/images/sphx_glr_plot_minimum_volume_level_sets_006.png
   :class: sphx-glr-single-img





.. GENERATED FROM PYTHON SOURCE LINES 176-180

Draw minimum volume level set in 2D
-----------------------------------

In this paragraph, we compute the minimum volume level set of a bivariate distribution.

.. GENERATED FROM PYTHON SOURCE LINES 183-184

Create a gaussian

.. GENERATED FROM PYTHON SOURCE LINES 184-199

.. code-block:: Python

    corr = ot.CorrelationMatrix(2)
    corr[0, 1] = 0.2
    copula = ot.NormalCopula(corr)
    x1 = ot.Normal(-1.0, 1)
    x2 = ot.Normal(2, 1)
    x_funk = ot.ComposedDistribution([x1, x2], copula)

    # Create a second gaussian
    x1 = ot.Normal(1.0, 1)
    x2 = ot.Normal(-2, 1)
    x_punk = ot.ComposedDistribution([x1, x2], copula)

    # Mix the distributions
    mixture = ot.Mixture([x_funk, x_punk], [0.5, 1.0])








.. GENERATED FROM PYTHON SOURCE LINES 200-203

.. code-block:: Python

    graph = mixture.drawPDF()
    view = viewer.View(graph)




.. image-sg:: /auto_probabilistic_modeling/distributions/images/sphx_glr_plot_minimum_volume_level_sets_007.png
   :alt: [X0,X1] iso-PDF
   :srcset: /auto_probabilistic_modeling/distributions/images/sphx_glr_plot_minimum_volume_level_sets_007.png
   :class: sphx-glr-single-img





.. GENERATED FROM PYTHON SOURCE LINES 204-205

For a multivariate distribution (with dimension greater than 1), the `computeMinimumVolumeLevelSetWithThreshold` uses Monte-Carlo sampling.

.. GENERATED FROM PYTHON SOURCE LINES 207-211

.. code-block:: Python

    ot.ResourceMap.SetAsUnsignedInteger(
        "Distribution-MinimumVolumeLevelSetSamplingSize", 1000
    )








.. GENERATED FROM PYTHON SOURCE LINES 212-215

We want to compute the minimum volume LevelSet which contains `alpha`=90% of the distribution.
The `threshold` is the value of the PDF corresponding the `alpha`-probability:
the points contained in the LevelSet have a PDF value lower or equal to this threshold.

.. GENERATED FROM PYTHON SOURCE LINES 217-222

.. code-block:: Python

    alpha = 0.9
    levelSet, threshold = mixture.computeMinimumVolumeLevelSetWithThreshold(alpha)
    threshold






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

 .. code-block:: none


    0.007697813348147185



.. GENERATED FROM PYTHON SOURCE LINES 223-257

.. code-block:: Python

    def drawLevelSetContour2D(
        distribution, numberOfPointsInXAxis, alpha, threshold, sampleSize=500
    ):
        """
        Compute the minimum volume LevelSet of measure equal to alpha and get the
        corresponding density value (named threshold).
        Generate a sample of the distribution and draw it.
        Draw a contour plot for the distribution, where the PDF is equal to threshold.
        """
        sample = distribution.getSample(sampleSize)
        X1min = sample[:, 0].getMin()[0]
        X1max = sample[:, 0].getMax()[0]
        X2min = sample[:, 1].getMin()[0]
        X2max = sample[:, 1].getMax()[0]
        xx = ot.Box([numberOfPointsInXAxis], ot.Interval([X1min], [X1max])).generate()
        yy = ot.Box([numberOfPointsInXAxis], ot.Interval([X2min], [X2max])).generate()
        xy = ot.Box(
            [numberOfPointsInXAxis, numberOfPointsInXAxis],
            ot.Interval([X1min, X2min], [X1max, X2max]),
        ).generate()
        data = distribution.computePDF(xy)
        graph = ot.Graph("", "X1", "X2", True, "upper right")
        labels = ["%.2f%%" % (100 * alpha)]
        contour = ot.Contour(xx, yy, data, [threshold], labels)
        contour.setColor("black")
        graph.setTitle(
            "%.2f%% of the distribution, sample size = %d" % (100 * alpha, sampleSize)
        )
        graph.add(contour)
        cloud = ot.Cloud(sample)
        graph.add(cloud)
        return graph









.. GENERATED FROM PYTHON SOURCE LINES 258-259

The following plot shows that 90% of the sample is contained in the LevelSet.

.. GENERATED FROM PYTHON SOURCE LINES 261-266

.. code-block:: Python

    numberOfPointsInXAxis = 50
    graph = drawLevelSetContour2D(mixture, numberOfPointsInXAxis, alpha, threshold)
    view = viewer.View(graph)
    plt.show()




.. image-sg:: /auto_probabilistic_modeling/distributions/images/sphx_glr_plot_minimum_volume_level_sets_008.png
   :alt: 90.00% of the distribution, sample size = 500
   :srcset: /auto_probabilistic_modeling/distributions/images/sphx_glr_plot_minimum_volume_level_sets_008.png
   :class: sphx-glr-single-img





.. GENERATED FROM PYTHON SOURCE LINES 267-268

Reset default settings

.. GENERATED FROM PYTHON SOURCE LINES 268-269

.. code-block:: Python

    ot.ResourceMap.Reload()








.. _sphx_glr_download_auto_probabilistic_modeling_distributions_plot_minimum_volume_level_sets.py:

.. only:: html

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

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

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

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

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