Note

Click here to download the full example code

Create unions or intersections of events¶

Abstract¶

This example illustrates system events, which are defined as unions or intersections of other events. We will show how to estimate their probability both with Monte-Carlo sampling (using the class ProbabilitySimulationAlgorithm) and with a first order approximation (using the class SystemFORM).

from __future__ import print_function
import openturns as ot
import openturns.viewer as otv
from matplotlib import pylab as plt
ot.Log.Show(ot.Log.NONE)

Intersection

The event defined as the intersection of several events is realized when all sub-events occurs:

$E_{sys} = \bigcap_{i=1}^N E_i$

Union

The event defined as the union of several events is realized when at least one sub-event occurs:

$E_{sys} = \bigcup_{i=1}^N E_i$

We consider a bivariate standard Gaussian random vector X = (X_1, X_2).

dim = 2
distribution = ot.Normal(dim)
X = ot.RandomVector(distribution)

We want to estimate the probability given by

$P = \mathbb{E}[\mathbf{1}_{\mathrm{Event}}(X_1, X_2)].$

We now build several events using intersections and unions.

We consider three functions f1, f2 and f3 :

f1 = ot.SymbolicFunction(['x0', 'x1'], ['x0'])
f2 = ot.SymbolicFunction(['x0', 'x1'], ['x1'])
f3 = ot.SymbolicFunction(['x0', 'x1'], ['x0+x1'])

We build CompositeRandomVector from these functions and the initial distribution.

Y1 = ot.CompositeRandomVector(f1, X)
Y2 = ot.CompositeRandomVector(f2, X)
Y3 = ot.CompositeRandomVector(f3, X)

We define three basic events $E_1=\{(x_0,x_1)~:~x_0 < 0 \}$ , $E_2=\{(x_0,x_1)~:~x_1 > 0 \}$ and $E_3=\{(x_0,x_1)~:~x_0+x_1>0 \}$ .

e1 = ot.ThresholdEvent(Y1, ot.Less(), 0.0)
e2 = ot.ThresholdEvent(Y2, ot.Greater(), 0.0)
e3 = ot.ThresholdEvent(Y3, ot.Greater(), 0.0)

The restriction of the domain $E_1$ to $[-4,4] \times [-4, 4]$ is the grey area.

myGraph = ot.Graph(r'Representation of the event $E_1$',
                   r'$x_1$', r'$x_2$', True, '')
data = [[-4, -4], [0, -4], [0, 4], [-4, 4]]
myPolygon = ot.Polygon(data)
myPolygon.setColor('grey')
myPolygon.setEdgeColor('black')
myGraph.add(myPolygon)
view = otv.View(myGraph)
axes = view.getAxes()
_ = axes[0].set_xlim(-4.0, 4.0)
_ = axes[0].set_ylim(-4.0, 4.0)

The restriction of the domain $E_2$ to $[-4,4] \times [-4, 4]$ is the grey area.

myGraph = ot.Graph(r'Representation of the event $E_2$',
                   r'$x_1$', r'$x_2$', True, '')
data = [[-4, 0], [4, 0], [4, 4], [-4, 4]]
myPolygon = ot.Polygon(data)
myPolygon.setColor('grey')
myPolygon.setEdgeColor('black')
myGraph.add(myPolygon)
view = otv.View(myGraph)
axes = view.getAxes()
_ = axes[0].set_xlim(-4.0, 4.0)
_ = axes[0].set_ylim(-4.0, 4.0)

The restriction of the domain $E_3$ to $[-4,4] \times [-4, 4]$ is the grey area.

myGraph = ot.Graph(r'Representation of the event $E_3$',
                   r'$x_1$', r'$x_2$', True, '')
data = [[-4, 4], [4, -4], [4, 4]]
myPolygon = ot.Polygon(data)
myPolygon.setColor('grey')
myPolygon.setEdgeColor('black')
myGraph.add(myPolygon)
view = otv.View(myGraph)
axes = view.getAxes()
_ = axes[0].set_xlim(-4.0, 4.0)
_ = axes[0].set_ylim(-4.0, 4.0)

We can define the intersection $E_4 = E_1 \bigcap E_2$ : that is the upper left quadrant.

e4 = ot.IntersectionEvent([e1, e2])

The restriction of the domain $E_4$ to $[-4,4] \times [-4, 4]$ is the grey area.

myGraph = ot.Graph(
    r'Representation of the event $E_4  = E_1 \bigcap E_2$', r'$x_1$', r'$x_2$', True, '')
data = [[-4, 0], [0, 0], [0, 4], [-4, 4]]
myPolygon = ot.Polygon(data)
myPolygon.setColor('grey')
myPolygon.setEdgeColor('black')
myGraph.add(myPolygon)
view = otv.View(myGraph)
axes = view.getAxes()
_ = axes[0].set_xlim(-4.0, 4.0)
_ = axes[0].set_ylim(-4.0, 4.0)

$Representation of the event $E_4 = E_1 \bigcap E_2$$

The probability of that event is $P_{E_4} = 1/4$ . A basic estimator is:

print("Probability of e4 : %.4f" % e4.getSample(10000).computeMean()[0])

Out:

Probability of e4 : 0.2504

We define the union $E_5 = E1 \bigcup E_2$ . It is the whole plan without the lower right quadrant.

e5 = ot.UnionEvent([e1, e2])

The restriction of the domain $E_5$ to $[-4,4] \times [-4, 4]$ is the grey area.

myGraph = ot.Graph(
    r'Representation of the event $E_5  = E_1 \bigcup E_2$', r'$x_1$', r'$x_2$', True, '')
data = [[-4, -4], [0, -4], [0, 0], [4, 0], [4, 4], [-4, 4]]
myPolygon = ot.Polygon(data)
myPolygon.setColor('grey')
myPolygon.setEdgeColor('black')
myGraph.add(myPolygon)
view = otv.View(myGraph)
axes = view.getAxes()
_ = axes[0].set_xlim(-4.0, 4.0)
_ = axes[0].set_ylim(-4.0, 4.0)

$Representation of the event $E_5 = E_1 \bigcup E_2$$

The probability of that event is $P_{E_5} = 3/4$ . A basic estimator is:

print("Probability of e5 : %.4f" % e5.getSample(10000).computeMean()[0])

Out:

Probability of e5 : 0.7489

It supports recursion. Let’s define $E_6 = E_1 \bigcup (E_2 \bigcap E_3)$ .

e6 = ot.UnionEvent([e1, ot.IntersectionEvent([e2, e3])])

First we draw the domain $E_6 = E_1 \bigcup (E_2 \bigcap E_3)$ :

myGraph = ot.Graph(
    r'Representation of the event $E_2 \bigcap E_3 $', r'$x_1$', r'$x_2$', True, '')
data = [[-4, 4], [0, 0], [4, 0], [4, 4]]
myPolygon = ot.Polygon(data)
myPolygon.setColor('grey')
myPolygon.setEdgeColor('black')
myGraph.add(myPolygon)
view = otv.View(myGraph)
axes = view.getAxes()
_ = axes[0].set_xlim(-4.0, 4.0)
_ = axes[0].set_ylim(-4.0, 4.0)

$Representation of the event $E_2 \bigcap E_3 $$

From the previous figures we easily deduce that the event $E_6 = E_1 \bigcup (E_2 \bigcap E_3)$ is the event $E_5$ and the probability is $P_{E_6} = 3/4$ . We can use a basic estimator and get :

print("Probability of e6 : %.4f" % e6.getSample(10000).computeMean()[0])

Out:

Probability of e6 : 0.7510

Usage with a Monte-Carlo algorithm¶

Of course, we can use simulation algorithms with this kind of events.

We set up a MonteCarloExperiment and a ProbabilitySimulationAlgorithm on the event $E_6$ .

experiment = ot.MonteCarloExperiment()
algo = ot.ProbabilitySimulationAlgorithm(e6, experiment)
algo.setMaximumOuterSampling(2500)
algo.setBlockSize(4)
algo.setMaximumCoefficientOfVariation(-1.0)
algo.run()

We retrieve the results and display the approximate probability and a confidence interval :

result = algo.getResult()
prb = result.getProbabilityEstimate()
print("Probability of e6 through MC : %.4f" % prb)
cl = result.getConfidenceLength()
print("Confidence interval MC : [%.4f, %.4f]" % (prb-0.5*cl, prb+0.5*cl))

Out:

Probability of e6 through MC : 0.7483
Confidence interval MC : [0.7388, 0.7578]

Usage with SystemFORM¶

The SystemFORM class implements an approximation method suitable for system events. The event must be in its disjunctive normal form (union of intersections, or a single intersection).

For system events, we always have to use the same root cause (input distribution). Here we use input variables with a normal distribution specified by its mean, standard deviation and correlation matrix.

dim = 5
mean = [200.0] * dim
mean[-1] = 60
mean[-2] = 60
sigma = [30.0] * dim
sigma[-1] = 15.0
R = ot.CorrelationMatrix(dim)
for i in range(dim):
    for j in range(i):
        R[i, j] = 0.5
dist = ot.Normal(mean, sigma, R)

As usual we create a RandomVector out of the input distribution.

X = ot.RandomVector(dist)

We define the leaf events thanks to SymbolicFunction.

inputs = ['M1', 'M2', 'M3', 'M4', 'M5']
e0 = ot.ThresholdEvent(ot.CompositeRandomVector(
    ot.SymbolicFunction(inputs, ['M1-M2+M4']), X), ot.Less(), 0.0)
e1 = ot.ThresholdEvent(ot.CompositeRandomVector(
    ot.SymbolicFunction(inputs, ['M2+2*M3-M4']), X), ot.Less(), 0.0)
e2 = ot.ThresholdEvent(ot.CompositeRandomVector(
    ot.SymbolicFunction(inputs, ['2*M3-2*M4-M5']), X), ot.Less(), 0.0)
e3 = ot.ThresholdEvent(ot.CompositeRandomVector(ot.SymbolicFunction(
    inputs, ['-(M1+M2+M4+M5-5*10.0)']), X), ot.Less(), 0.0)
e4 = ot.ThresholdEvent(ot.CompositeRandomVector(
    ot.SymbolicFunction(inputs, ['-(M2+2*M3+M4-5*40.0)']), X), ot.Less(), 0.0)

We consider a system event in disjunctive normal form (union of intersections):

event = ot.UnionEvent([ot.IntersectionEvent([e0, e3, e4]),
                       ot.IntersectionEvent([e2, e3, e4])])

We can estimate the probability of the event with basic sampling.

print("Probability of the event : %.4f" %
      event.getSample(10000).computeMean()[0])

Out:

Probability of the event : 0.0788

We can also run a systemFORM algorithm to estimate the probability differently.

We first set up a solver to find the design point.

solver = ot.AbdoRackwitz()
solver.setMaximumIterationNumber(1000)
solver.setMaximumAbsoluteError(1.0e-3)
solver.setMaximumRelativeError(1.0e-3)
solver.setMaximumResidualError(1.0e-3)
solver.setMaximumConstraintError(1.0e-3)

We build the SystemFORM algorithm from the solver, the event and a starting point (here the mean) and then run the algorithm.

algo = ot.SystemFORM(solver, event, mean)
algo.run()

We store the result and display the probability.

result = algo.getResult()
prbSystemFORM = result.getEventProbability()
print("Probability of the event (SystemFORM) : %.4f" % prbSystemFORM)

Out:

Probability of the event (SystemFORM) : 0.0788

Display all figures

plt.show()

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

Gallery generated by Sphinx-Gallery

OpenTURNS

An Open source initiative for the Treatment of Uncertainties, Risks'N Statistics

Table of Contents

Previous topic

Next topic

This Page

Create unions or intersections of events¶

Abstract¶

Usage with a Monte-Carlo algorithm¶

Usage with SystemFORM¶