Compare unconditional and conditional histogramsΒΆ

In this example, we compare unconditional and conditional histograms for a simulation. We consider the flooding model. Let g be a function which takes four inputs Q, K_s, Z_v and Z_m and returns one output S.

We first consider the (unconditional) distribution of the input Q.

Let t be a given threshold on the output S: we consider the event S > t. Then we consider the conditional distribution of the input Q given that S > t that is to say Q|S > t.

If these two distributions are significantly different, we conclude that the input Q has an impact on the event S > t.

In order to approximate the distribution of the output S, we perform a Monte-Carlo simulation with size 500. The threshold t is chosen as the 90% quantile of the empirical distribution of S. In this example, the distribution is aproximated by its empirical histogram (but this could be done with another distribution approximation as well, such as kernel smoothing for example).

import numpy as np
from openturns.usecases import flood_model
import openturns as ot
import openturns.viewer as viewer
from matplotlib import pylab as plt

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

We use the FloodModel data class that contains all the case parameters.

fm = flood_model.FloodModel()

Create an input sample from the joint distribution defined in the data class. We build an output sample by taking the image by the model.

size = 500
inputSample = fm.distribution.getSample(size)
inputSample[:5]
Q (m3/s)KsZv (m)Zm (m)B (m)L (m)Zb (m)Hd (m)
02032.97828.1643149.8182354.44882298.09834997.51155.276753.987806
1831.178432.0659849.857854.29531298.31574997.29755.187412.030507
21741.77619.3668149.0897555.0745299.14334999.43255.666932.719918
3800.47640.0074349.1621655.03673299.59985002.71255.557153.080748
4917.983538.2301849.1987854.97124302.27655008.60755.366593.816204


outputSample = fm.model(inputSample)
outputSample[:5]
HSC
03.470401-5.9759261.034781
11.900478-5.4596431.140406
23.659279-5.6378221.102684
31.492342-7.9833980.7745734
41.66541-8.31860.7507753


Merge the input and output samples into a single sample.

sample = ot.Sample(inputSample)
sample.stack(outputSample)
sample[0:5]
Q (m3/s)KsZv (m)Zm (m)B (m)L (m)Zb (m)Hd (m)HSC
02032.97828.1643149.8182354.44882298.09834997.51155.276753.9878063.470401-5.9759261.034781
1831.178432.0659849.857854.29531298.31574997.29755.187412.0305071.900478-5.4596431.140406
21741.77619.3668149.0897555.0745299.14334999.43255.666932.7199183.659279-5.6378221.102684
3800.47640.0074349.1621655.03673299.59985002.71255.557153.0807481.492342-7.9833980.7745734
4917.983538.2301849.1987854.97124302.27655008.60755.366593.8162041.66541-8.31860.7507753


Extract the first column of inputSample into the sample of the flowrates Q.

sampleQ = inputSample[:, 0]

The next cell defines a function that computes the conditional sample of a component given that the a marginal (defined by its index criteriaComponent) exceeds a given threshold, defined by its quantile level.

def computeConditionnedSample(
    sample, alpha=0.9, criteriaComponent=None, selectedComponent=0
):
    """
    Return values from the selectedComponent-th component of the sample.
    Selects the values according to the alpha-level quantile of
    the criteriaComponent-th component of the sample.
    """
    dim = sample.getDimension()
    if criteriaComponent is None:
        criteriaComponent = dim - 1
    sortedSample = sample.sortAccordingToAComponent(criteriaComponent)
    quantiles = sortedSample.computeQuantilePerComponent(alpha)
    quantileValue = quantiles[criteriaComponent]
    sortedSampleCriteria = sortedSample[:, criteriaComponent]
    indices = np.where(np.array(sortedSampleCriteria.asPoint()) > quantileValue)[0]
    conditionnedSortedSample = sortedSample[int(indices[0]):, selectedComponent]
    return conditionnedSortedSample

Create an histogram for the unconditional flowrates.

numberOfBins = 10
histogram = ot.HistogramFactory().buildAsHistogram(sampleQ, numberOfBins)

Extract the sub-sample of the input flowrates Q which leads to large values of the output S.

Search the index of the marginal S in the columns of the sample.

criteriaComponent = list(sample.getDescription()).index("S")
criteriaComponent
9
alpha = 0.9
selectedComponent = 0
conditionnedSampleQ = computeConditionnedSample(
    sample, alpha, criteriaComponent, selectedComponent
)

We could as well use:

# conditionnedHistogram = ot.HistogramFactory().buildAsHistogram(conditionnedSampleQ)

but this creates an histogram with new classes, corresponding to conditionnedSampleQ. We want to use exactly the same classes as the full sample, so that the two histograms match.

first = histogram.getFirst()
width = histogram.getWidth()
conditionnedHistogram = ot.HistogramFactory().buildAsHistogram(
    conditionnedSampleQ, first, width
)

Then creates a graphics with the unconditional and the conditional histograms.

graph = histogram.drawPDF()
graph.setLegends(["Q"])
#
graphConditionnalQ = conditionnedHistogram.drawPDF()
graphConditionnalQ.setColors(["blue"])
graphConditionnalQ.setLegends([r"$Q | S > S_{%s}$" % (alpha)])
graph.add(graphConditionnalQ)
view = viewer.View(graph)

plt.show()
Q (m3/s) PDF

We see that the two histograms are very different. The high values of the input Q seem to often lead to a high value of the output S.

We could explore this situation further by comparing the unconditional distribution of Q (which is known in this case) with the conditonal distribution of Q | S > t, estimated by kernel smoothing. This would have the advantage of accuracy, since the kernel smoothing is a more accurate approximation of a distribution than the histogram.