Draw minimum volume level sets

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

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

Draw minimum volume level set in 1D

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

With a Normal, minimum volume LevelSet

n = ot.Normal()
graph = n.drawPDF()
view = viewer.View(graph)
plot minimum volume level sets

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.

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

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

levelSet.contains([0.0])
True
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
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
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
graph = drawLevelSet1D(n, levelSet, alpha, threshold)
view = viewer.View(graph)
90.00% of the distribution, sample size = 100,

With a Normal, minimum volume Interval

interval = n.computeMinimumVolumeInterval(alpha)
interval

[-1.64485, 1.64485]



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

The computeMinimumVolumeInterval returns an Interval.

graph = drawPDFAndInterval1D(n, interval, alpha)
view = viewer.View(graph)
90.00% of the distribution, lower bound = -1.645, upper bound = 1.645

With a Mixture, minimum volume LevelSet

m = ot.Mixture([ot.Normal(-5.0, 1.0), ot.Normal(5.0, 1.0)], [0.2, 0.8])
graph = m.drawPDF()
view = viewer.View(graph)
plot minimum volume level sets
alpha = 0.9
levelSet, threshold = m.computeMinimumVolumeLevelSetWithThreshold(alpha)
threshold
0.04667473141153258

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

graph = drawLevelSet1D(m, levelSet, alpha, threshold, 1000)
view = viewer.View(graph)
90.00% of the distribution, sample size = 1000,

With a Mixture, minimum volume Interval

interval = m.computeMinimumVolumeInterval(alpha)
interval

[-5.44003, 6.72227]



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

graph = drawPDFAndInterval1D(m, interval, alpha)
view = viewer.View(graph)
90.00% of the distribution, lower bound = -5.440, upper bound = 6.722

Draw minimum volume level set in 2D

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

Create a gaussian

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])
graph = mixture.drawPDF()
view = viewer.View(graph)
[X0,X1] iso-PDF

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

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

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.

alpha = 0.9
levelSet, threshold = mixture.computeMinimumVolumeLevelSetWithThreshold(alpha)
threshold
0.0076265926873795705
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

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

numberOfPointsInXAxis = 50
graph = drawLevelSetContour2D(mixture, numberOfPointsInXAxis, alpha, threshold)
view = viewer.View(graph)
plt.show()
90.00% of the distribution, sample size = 500
plt.show()

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