.. 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-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 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-93 .. 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 94-97 .. code-block:: default 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 98-100 With a Normal, minimum volume Interval ^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^ .. GENERATED FROM PYTHON SOURCE LINES 102-106 .. code-block:: default interval = n.computeMinimumVolumeInterval(alpha) interval .. raw:: html

[-1.64485, 1.64485]



.. GENERATED FROM PYTHON SOURCE LINES 107-123 .. 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 124-125 The `computeMinimumVolumeInterval` returns an `Interval`. .. GENERATED FROM PYTHON SOURCE LINES 127-130 .. code-block:: default 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 131-133 With a Mixture, minimum volume LevelSet ^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^ .. GENERATED FROM PYTHON SOURCE LINES 135-137 .. code-block:: default m = ot.Mixture([ot.Normal(-5., 1.), ot.Normal(5., 1.)], [0.2, 0.8]) .. GENERATED FROM PYTHON SOURCE LINES 138-141 .. code-block:: default 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 142-146 .. 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 147-148 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 150-153 .. code-block:: default 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 154-156 With a Mixture, minimum volume Interval ^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^ .. GENERATED FROM PYTHON SOURCE LINES 158-161 .. code-block:: default interval = m.computeMinimumVolumeInterval(alpha) interval .. raw:: html

[-5.44003, 6.72227]



.. GENERATED FROM PYTHON SOURCE LINES 162-163 The `computeMinimumVolumeInterval` returns an `Interval`. The bounds of this interval are different from the previous `LevelSet`. .. GENERATED FROM PYTHON SOURCE LINES 165-169 .. code-block:: default 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 170-174 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 177-178 Create a gaussian .. GENERATED FROM PYTHON SOURCE LINES 178-193 .. 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 194-197 .. code-block:: default 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 198-199 For a multivariate distribution (with dimension greater than 1), the `computeMinimumVolumeLevelSetWithThreshold` uses Monte-Carlo sampling. .. GENERATED FROM PYTHON SOURCE LINES 201-204 .. code-block:: default ot.ResourceMap.SetAsUnsignedInteger( "Distribution-MinimumVolumeLevelSetSamplingSize", 1000) .. GENERATED FROM PYTHON SOURCE LINES 205-206 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 208-213 .. 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 214-245 .. 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 246-247 The following plot shows that 90% of the sample is contained in the LevelSet. .. GENERATED FROM PYTHON SOURCE LINES 249-253 .. code-block:: default 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 254-256 .. code-block:: default plt.show() .. rst-class:: sphx-glr-timing **Total running time of the script:** ( 0 minutes 0.621 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 `_