.. 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_probabilistic_modeling_distributions_plot_minimum_volume_level_set_1D.py:
Draw minimum volume level set in 1D
===================================
In this example, we compute the minimum volume level set of a univariate distribution.
.. code-block:: default
import openturns as ot
import openturns.viewer as viewer
from matplotlib import pylab as plt
ot.Log.Show(ot.Log.NONE)
With a Normal, minimum volume LevelSet
--------------------------------------
.. code-block:: default
n = ot.Normal()
.. code-block:: default
graph = n.drawPDF()
view = viewer.View(graph)
.. image:: /auto_probabilistic_modeling/distributions/images/sphx_glr_plot_minimum_volume_level_set_1D_001.png
:alt: plot minimum volume level set 1D
:class: sphx-glr-single-img
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.
.. code-block:: default
alpha = 0.9
levelSet, threshold = n.computeMinimumVolumeLevelSetWithThreshold(alpha)
threshold
.. rst-class:: sphx-glr-script-out
Out:
.. code-block:: none
0.1031356403794346
The `LevelSet` has a `contains` method. Obviously, the point 0 is in the LevelSet.
.. code-block:: default
levelSet.contains([0.])
.. rst-class:: sphx-glr-script-out
Out:
.. code-block:: none
True
.. code-block:: default
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
.. code-block:: default
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
.. code-block:: default
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
.. code-block:: default
graph = drawLevelSet1D(n, levelSet, alpha, threshold)
view = viewer.View(graph)
.. image:: /auto_probabilistic_modeling/distributions/images/sphx_glr_plot_minimum_volume_level_set_1D_002.png
:alt: 90.00% of the distribution, sample size = 100,
:class: sphx-glr-single-img
With a Normal, minimum volume Interval
--------------------------------------
.. code-block:: default
interval = n.computeMinimumVolumeInterval(alpha)
interval
.. raw:: html
[-1.64485, 1.64485]
.. code-block:: default
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.],[xmin,yvalue],[xmax,yvalue],[xmax,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
The `computeMinimumVolumeInterval` returns an `Interval`.
.. code-block:: default
graph = drawPDFAndInterval1D(n, interval, alpha)
view = viewer.View(graph)
.. image:: /auto_probabilistic_modeling/distributions/images/sphx_glr_plot_minimum_volume_level_set_1D_003.png
:alt: 90.00% of the distribution, lower bound = -1.645, upper bound = 1.645
:class: sphx-glr-single-img
With a Mixture, minimum volume LevelSet
---------------------------------------
.. code-block:: default
m = ot.Mixture([ot.Normal(-5.,1.),ot.Normal(5.,1.)],[0.2,0.8])
.. code-block:: default
graph = m.drawPDF()
view = viewer.View(graph)
.. image:: /auto_probabilistic_modeling/distributions/images/sphx_glr_plot_minimum_volume_level_set_1D_004.png
:alt: plot minimum volume level set 1D
:class: sphx-glr-single-img
.. code-block:: default
alpha = 0.9
levelSet, threshold = m.computeMinimumVolumeLevelSetWithThreshold(alpha)
threshold
.. rst-class:: sphx-glr-script-out
Out:
.. code-block:: none
0.04667473141178892
The interesting point is that a `LevelSet` may be non-contiguous. In the current mixture example, this is not an interval.
.. code-block:: default
graph = drawLevelSet1D(m, levelSet, alpha, threshold, 1000)
view = viewer.View(graph)
.. image:: /auto_probabilistic_modeling/distributions/images/sphx_glr_plot_minimum_volume_level_set_1D_005.png
:alt: 90.00% of the distribution, sample size = 1000,
:class: sphx-glr-single-img
With a Mixture, minimum volume Interval
---------------------------------------
.. code-block:: default
interval = m.computeMinimumVolumeInterval(alpha)
interval
.. raw:: html
[-5.44003, 6.72227]
The `computeMinimumVolumeInterval` returns an `Interval`. The bounds of this interval are different from the previous `LevelSet`.
.. code-block:: default
graph = drawPDFAndInterval1D(m, interval, alpha)
view = viewer.View(graph)
plt.show()
.. image:: /auto_probabilistic_modeling/distributions/images/sphx_glr_plot_minimum_volume_level_set_1D_006.png
:alt: 90.00% of the distribution, lower bound = -5.440, upper bound = 6.722
:class: sphx-glr-single-img
.. rst-class:: sphx-glr-timing
**Total running time of the script:** ( 0 minutes 0.502 seconds)
.. _sphx_glr_download_auto_probabilistic_modeling_distributions_plot_minimum_volume_level_set_1D.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_minimum_volume_level_set_1D.py `
.. container:: sphx-glr-download sphx-glr-download-jupyter
:download:`Download Jupyter notebook: plot_minimum_volume_level_set_1D.ipynb `
.. only:: html
.. rst-class:: sphx-glr-signature
`Gallery generated by Sphinx-Gallery `_