.. 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 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_sets.py: Draw minimum volume level sets ============================== .. GENERATED FROM PYTHON SOURCE LINES 6-11 .. code-block:: default 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 12-16 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 19-21 With a Normal, minimum volume LevelSet ^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^ .. GENERATED FROM PYTHON SOURCE LINES 23-25 .. code-block:: default n = ot.Normal() .. GENERATED FROM PYTHON SOURCE LINES 26-29 .. code-block:: default graph = n.drawPDF() view = viewer.View(graph) .. image:: /auto_probabilistic_modeling/distributions/images/sphx_glr_plot_minimum_volume_level_sets_001.png :alt: plot minimum volume level sets :class: sphx-glr-single-img .. GENERATED FROM PYTHON SOURCE LINES 30-31 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 33-37 .. code-block:: default alpha = 0.9 levelSet, threshold = n.computeMinimumVolumeLevelSetWithThreshold(alpha) threshold .. rst-class:: sphx-glr-script-out Out: .. code-block:: none 0.10313564037537128 .. GENERATED FROM PYTHON SOURCE LINES 38-39 The `LevelSet` has a `contains` method. Obviously, the point 0 is in the LevelSet. .. GENERATED FROM PYTHON SOURCE LINES 41-44 .. code-block:: default levelSet.contains([0.]) .. rst-class:: sphx-glr-script-out Out: .. code-block:: none True .. GENERATED FROM PYTHON SOURCE LINES 45-65 .. 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 .. GENERATED FROM PYTHON SOURCE LINES 66-77 .. 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 .. GENERATED FROM PYTHON SOURCE LINES 78-92 .. 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 .. GENERATED FROM PYTHON SOURCE LINES 93-96 .. 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_sets_002.png :alt: 90.00% of the distribution, sample size = 100, :class: sphx-glr-single-img .. GENERATED FROM PYTHON SOURCE LINES 97-99 With a Normal, minimum volume Interval ^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^ .. GENERATED FROM PYTHON SOURCE LINES 101-105 .. code-block:: default interval = n.computeMinimumVolumeInterval(alpha) interval .. raw:: html

[-1.64485, 1.64485]



.. GENERATED FROM PYTHON SOURCE LINES 106-121 .. 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 .. GENERATED FROM PYTHON SOURCE LINES 122-123 The `computeMinimumVolumeInterval` returns an `Interval`. .. GENERATED FROM PYTHON SOURCE LINES 125-128 .. code-block:: default graph = drawPDFAndInterval1D(n, interval, alpha) view = viewer.View(graph) .. image:: /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 :class: sphx-glr-single-img .. GENERATED FROM PYTHON SOURCE LINES 129-131 With a Mixture, minimum volume LevelSet ^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^ .. GENERATED FROM PYTHON SOURCE LINES 133-135 .. code-block:: default m = ot.Mixture([ot.Normal(-5.,1.),ot.Normal(5.,1.)],[0.2,0.8]) .. GENERATED FROM PYTHON SOURCE LINES 136-139 .. code-block:: default graph = m.drawPDF() view = viewer.View(graph) .. image:: /auto_probabilistic_modeling/distributions/images/sphx_glr_plot_minimum_volume_level_sets_004.png :alt: plot minimum volume level sets :class: sphx-glr-single-img .. GENERATED FROM PYTHON SOURCE LINES 140-144 .. code-block:: default alpha = 0.9 levelSet, threshold = m.computeMinimumVolumeLevelSetWithThreshold(alpha) threshold .. rst-class:: sphx-glr-script-out Out: .. code-block:: none 0.04667473141153258 .. GENERATED FROM PYTHON SOURCE LINES 145-146 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 148-151 .. 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_sets_005.png :alt: 90.00% of the distribution, sample size = 1000, :class: sphx-glr-single-img .. GENERATED FROM PYTHON SOURCE LINES 152-154 With a Mixture, minimum volume Interval ^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^ .. GENERATED FROM PYTHON SOURCE LINES 156-159 .. code-block:: default interval = m.computeMinimumVolumeInterval(alpha) interval .. raw:: html

[-5.44003, 6.72227]



.. GENERATED FROM PYTHON SOURCE LINES 160-161 The `computeMinimumVolumeInterval` returns an `Interval`. The bounds of this interval are different from the previous `LevelSet`. .. GENERATED FROM PYTHON SOURCE LINES 163-167 .. code-block:: default graph = drawPDFAndInterval1D(m, interval, alpha) view = viewer.View(graph) .. image:: /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 :class: sphx-glr-single-img .. GENERATED FROM PYTHON SOURCE LINES 168-172 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 175-176 Create a gaussian .. GENERATED FROM PYTHON SOURCE LINES 176-191 .. code-block:: default corr = ot.CorrelationMatrix(2) corr[0, 1] = 0.2 copula = ot.NormalCopula(corr) x1 = ot.Normal(-1., 1) x2 = ot.Normal(2, 1) x_funk = ot.ComposedDistribution([x1, x2], copula) # Create a second gaussian x1 = ot.Normal(1.,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.]) .. GENERATED FROM PYTHON SOURCE LINES 192-195 .. code-block:: default graph = mixture.drawPDF() view = viewer.View(graph) .. image:: /auto_probabilistic_modeling/distributions/images/sphx_glr_plot_minimum_volume_level_sets_007.png :alt: [X0,X1] iso-PDF :class: sphx-glr-single-img .. GENERATED FROM PYTHON SOURCE LINES 196-197 For a multivariate distribution (with dimension greater than 1), the `computeMinimumVolumeLevelSetWithThreshold` uses Monte-Carlo sampling. .. GENERATED FROM PYTHON SOURCE LINES 199-201 .. code-block:: default ot.ResourceMap.SetAsUnsignedInteger("Distribution-MinimumVolumeLevelSetSamplingSize",1000) .. GENERATED FROM PYTHON SOURCE LINES 202-203 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 205-210 .. code-block:: default alpha = 0.9 levelSet, threshold = mixture.computeMinimumVolumeLevelSetWithThreshold(alpha) threshold .. rst-class:: sphx-glr-script-out Out: .. code-block:: none 0.0076863340815168865 .. GENERATED FROM PYTHON SOURCE LINES 211-241 .. code-block:: default 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, 'topright') 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 242-243 The following plot shows that 90% of the sample is contained in the LevelSet. .. GENERATED FROM PYTHON SOURCE LINES 245-249 .. code-block:: default numberOfPointsInXAxis = 50 graph = drawLevelSetContour2D(mixture, numberOfPointsInXAxis, alpha, threshold) view = viewer.View(graph) plt.show() .. image:: /auto_probabilistic_modeling/distributions/images/sphx_glr_plot_minimum_volume_level_sets_008.png :alt: 90.00% of the distribution, sample size = 500 :class: sphx-glr-single-img .. GENERATED FROM PYTHON SOURCE LINES 250-252 .. code-block:: default plt.show() .. rst-class:: sphx-glr-timing **Total running time of the script:** ( 0 minutes 1.008 seconds) .. _sphx_glr_download_auto_probabilistic_modeling_distributions_plot_minimum_volume_level_sets.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_sets.py ` .. container:: sphx-glr-download sphx-glr-download-jupyter :download:`Download Jupyter notebook: plot_minimum_volume_level_sets.ipynb ` .. only:: html .. rst-class:: sphx-glr-signature `Gallery generated by Sphinx-Gallery `_