Estimate threshold exceedance iteratively

In this example, we use the IterativeThresholdExceedance class to count the number of threshold exceedances.

import openturns as ot
import openturns.viewer as otv

We first create a one-dimensional Gaussian random variable to generate data.

dim = 1
distNormal = ot.Normal(dim)

Let us consider a threshold value of 1.0. Each data value higher than 1.0 is counted as one exceedance. The counter used by the IterativeThresholdExceedance class is updated iteratively.

thresholdValue = 1.0
iterThreshold = ot.IterativeThresholdExceedance(dim, thresholdValue)

A simple computation shows that the probability of the data being higher than is (with 4 significant digits).

distribution = ot.Normal()
exactProbability = distribution.computeComplementaryCDF(thresholdValue)
print("Exact probability: ", exactProbability)

Out:

Exact probability:  0.15865525393145702

We can now perform the simulations. In our case most of the data fall below the specified threshold value so the number of exceedances should be low.

We first increment the object one Point at a time. At any given step the current number of exceedance is obtained with the getThresholdExceedance() method.

size = 5000
exceedanceNumbers = ot.Sample()
probabilityEstimateSample = ot.Sample()
for i in range(size):
    point = distNormal.getRealization()
    iterThreshold.increment(point)
    numberOfExceedances = iterThreshold.getThresholdExceedance()[0]
    exceedanceNumbers.add([numberOfExceedances])
    probabilityEstimate = numberOfExceedances / iterThreshold.getIterationNumber()
    probabilityEstimateSample.add([probabilityEstimate])

We display the evolution of the number of exceedances.

curve = ot.Curve(exceedanceNumbers)
curve.setLegend("number of exceedance")
#
graph = ot.Graph(
    "Evolution of the number of exceedances",
    "iteration nb",
    "number of exceedances",
    True,
)
graph.add(curve)
graph.setLegends(["number of exceedances"])
graph.setLegendPosition("bottomright")
view = otv.View(graph)

The following plot shows that the probability of exceeding the threshold converges.

palette = ot.Drawable().BuildDefaultPalette(2)
iterationSample = ot.Sample.BuildFromPoint(range(1, size + 1))
curve = ot.Curve(iterationSample, probabilityEstimateSample)
curve.setLegend("Prob. of exceeding the threshold")
curve.setColor(palette[0])
#
exactCurve = ot.Curve([1, size], [exactProbability, exactProbability])
exactCurve.setLegend("Exact")
exactCurve.setColor(palette[1])
#
graph = ot.Graph(
    "Evolution of the sample probability",
    "iteration nb",
    "estimate of the probability",
    True,
)
graph.add(curve)
graph.add(exactCurve)
graph.setLegendPosition("topleft")
graph.setLogScale(ot.GraphImplementation.LOGX)
view = otv.View(graph)

We can also increment with a Sample.

sample = distNormal.getSample(size)
iterThreshold.increment(sample)
numberOfExceedances = iterThreshold.getThresholdExceedance()[0]
print("Number of exceedance: ", numberOfExceedances)

# The empirical probability is close to the exact value.
numberOfExceedances = iterThreshold.getThresholdExceedance()[0]
probabilityEstimate = numberOfExceedances / iterThreshold.getIterationNumber()
print("Empirical exceedance prb: ", probabilityEstimate)

otv.View.ShowAll()

Out:

Number of exceedance:  1556.0
Empirical exceedance prb:  0.1556