.. DO NOT EDIT. .. THIS FILE WAS AUTOMATICALLY GENERATED BY SPHINX-GALLERY. .. TO MAKE CHANGES, EDIT THE SOURCE PYTHON FILE: .. "auto_reliability_sensitivity/reliability/plot_event_system.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_reliability_sensitivity_reliability_plot_event_system.py: Create unions or intersections of events ======================================== .. GENERATED FROM PYTHON SOURCE LINES 6-11 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 :class:`~openturns.ProbabilitySimulationAlgorithm`) and with a first order approximation (using the class :class:`~openturns.SystemFORM`). .. GENERATED FROM PYTHON SOURCE LINES 11-18 .. code-block:: default 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) .. GENERATED FROM PYTHON SOURCE LINES 19-33 **Intersection** The event defined as the intersection of several events is realized when all sub-events occurs: .. math:: 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: .. math:: E_{sys} = \bigcup_{i=1}^N E_i .. GENERATED FROM PYTHON SOURCE LINES 36-37 We consider a bivariate standard Gaussian random vector `X = (X_1, X_2)`. .. GENERATED FROM PYTHON SOURCE LINES 37-42 .. code-block:: default dim = 2 distribution = ot.Normal(dim) X = ot.RandomVector(distribution) .. GENERATED FROM PYTHON SOURCE LINES 43-51 We want to estimate the probability given by .. math:: P = \mathbb{E}[\mathbf{1}_{\mathrm{Event}}(X_1, X_2)]. We now build several events using intersections and unions. .. GENERATED FROM PYTHON SOURCE LINES 53-54 We consider three functions `f1`, `f2` and `f3` : .. GENERATED FROM PYTHON SOURCE LINES 54-58 .. code-block:: default f1 = ot.SymbolicFunction(['x0', 'x1'], ['x0']) f2 = ot.SymbolicFunction(['x0', 'x1'], ['x1']) f3 = ot.SymbolicFunction(['x0', 'x1'], ['x0+x1']) .. GENERATED FROM PYTHON SOURCE LINES 59-60 We build :class:`~openturns.CompositeRandomVector` from these functions and the initial distribution. .. GENERATED FROM PYTHON SOURCE LINES 60-64 .. code-block:: default Y1 = ot.CompositeRandomVector(f1, X) Y2 = ot.CompositeRandomVector(f2, X) Y3 = ot.CompositeRandomVector(f3, X) .. GENERATED FROM PYTHON SOURCE LINES 65-66 We define three basic events :math:`E_1=\{(x_0,x_1)~:~x_0 < 0 \}`, :math:`E_2=\{(x_0,x_1)~:~x_1 > 0 \}` and :math:`E_3=\{(x_0,x_1)~:~x_0+x_1>0 \}`. .. GENERATED FROM PYTHON SOURCE LINES 66-70 .. code-block:: default e1 = ot.ThresholdEvent(Y1, ot.Less(), 0.0) e2 = ot.ThresholdEvent(Y2, ot.Greater(), 0.0) e3 = ot.ThresholdEvent(Y3, ot.Greater(), 0.0) .. GENERATED FROM PYTHON SOURCE LINES 71-72 The restriction of the domain :math:`E_1` to :math:`[-4,4] \times [-4, 4]` is the grey area. .. GENERATED FROM PYTHON SOURCE LINES 72-84 .. code-block:: default 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) .. image:: /auto_reliability_sensitivity/reliability/images/sphx_glr_plot_event_system_001.png :alt: Representation of the event $E_1$ :class: sphx-glr-single-img .. GENERATED FROM PYTHON SOURCE LINES 85-86 The restriction of the domain :math:`E_2` to :math:`[-4,4] \times [-4, 4]` is the grey area. .. GENERATED FROM PYTHON SOURCE LINES 86-97 .. code-block:: default 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) .. image:: /auto_reliability_sensitivity/reliability/images/sphx_glr_plot_event_system_002.png :alt: Representation of the event $E_2$ :class: sphx-glr-single-img .. GENERATED FROM PYTHON SOURCE LINES 98-99 The restriction of the domain :math:`E_3` to :math:`[-4,4] \times [-4, 4]` is the grey area. .. GENERATED FROM PYTHON SOURCE LINES 99-110 .. code-block:: default 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) .. image:: /auto_reliability_sensitivity/reliability/images/sphx_glr_plot_event_system_003.png :alt: Representation of the event $E_3$ :class: sphx-glr-single-img .. GENERATED FROM PYTHON SOURCE LINES 111-112 We can define the intersection :math:`E_4 = E_1 \bigcap E_2`: that is the upper left quadrant. .. GENERATED FROM PYTHON SOURCE LINES 112-114 .. code-block:: default e4 = ot.IntersectionEvent([e1, e2]) .. GENERATED FROM PYTHON SOURCE LINES 115-116 The restriction of the domain :math:`E_4` to :math:`[-4,4] \times [-4, 4]` is the grey area. .. GENERATED FROM PYTHON SOURCE LINES 116-127 .. code-block:: default 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) .. image:: /auto_reliability_sensitivity/reliability/images/sphx_glr_plot_event_system_004.png :alt: Representation of the event $E_4 = E_1 \bigcap E_2$ :class: sphx-glr-single-img .. GENERATED FROM PYTHON SOURCE LINES 128-129 The probability of that event is :math:`P_{E_4} = 1/4`. A basic estimator is: .. GENERATED FROM PYTHON SOURCE LINES 129-131 .. code-block:: default print("Probability of e4 : %.4f"%e4.getSample(10000).computeMean()[0] ) .. rst-class:: sphx-glr-script-out Out: .. code-block:: none Probability of e4 : 0.2523 .. GENERATED FROM PYTHON SOURCE LINES 132-133 We define the union :math:`E_5 = E1 \bigcup E_2`. It is the whole plan without the lower right quadrant. .. GENERATED FROM PYTHON SOURCE LINES 133-135 .. code-block:: default e5 = ot.UnionEvent([e1, e2]) .. GENERATED FROM PYTHON SOURCE LINES 136-137 The restriction of the domain :math:`E_5` to :math:`[-4,4] \times [-4, 4]` is the grey area. .. GENERATED FROM PYTHON SOURCE LINES 137-148 .. code-block:: default 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) .. image:: /auto_reliability_sensitivity/reliability/images/sphx_glr_plot_event_system_005.png :alt: Representation of the event $E_5 = E_1 \bigcup E_2$ :class: sphx-glr-single-img .. GENERATED FROM PYTHON SOURCE LINES 149-150 The probability of that event is :math:`P_{E_5} = 3/4`. A basic estimator is: .. GENERATED FROM PYTHON SOURCE LINES 150-152 .. code-block:: default print("Probability of e5 : %.4f"%e5.getSample(10000).computeMean()[0] ) .. rst-class:: sphx-glr-script-out Out: .. code-block:: none Probability of e5 : 0.7498 .. GENERATED FROM PYTHON SOURCE LINES 153-154 It supports recursion. Let's define :math:`E_6 = E_1 \bigcup (E_2 \bigcap E_3)`. .. GENERATED FROM PYTHON SOURCE LINES 154-157 .. code-block:: default e6 = ot.UnionEvent([e1, ot.IntersectionEvent([e2, e3])]) .. GENERATED FROM PYTHON SOURCE LINES 158-159 First we draw the domain :math:`E_6 = E_1 \bigcup (E_2 \bigcap E_3)` : .. GENERATED FROM PYTHON SOURCE LINES 159-170 .. code-block:: default 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) .. image:: /auto_reliability_sensitivity/reliability/images/sphx_glr_plot_event_system_006.png :alt: Representation of the event $E_2 \bigcap E_3 $ :class: sphx-glr-single-img .. GENERATED FROM PYTHON SOURCE LINES 171-173 From the previous figures we easily deduce that the event :math:`E_6 = E_1 \bigcup (E_2 \bigcap E_3)` is the event :math:`E_5` and the probability is :math:`P_{E_6} = 3/4`. We can use a basic estimator and get : .. GENERATED FROM PYTHON SOURCE LINES 173-176 .. code-block:: default print("Probability of e6 : %.4f"%e6.getSample(10000).computeMean()[0] ) .. rst-class:: sphx-glr-script-out Out: .. code-block:: none Probability of e6 : 0.7487 .. GENERATED FROM PYTHON SOURCE LINES 177-181 Usage with a Monte-Carlo algorithm ---------------------------------- Of course, we can use simulation algorithms with this kind of events. .. GENERATED FROM PYTHON SOURCE LINES 183-184 We set up a :class:`~openturns.MonteCarloExperiment` and a :class:`~openturns.ProbabilitySimulationAlgorithm` on the event :math:`E_6`. .. GENERATED FROM PYTHON SOURCE LINES 184-191 .. code-block:: default experiment = ot.MonteCarloExperiment() algo = ot.ProbabilitySimulationAlgorithm(e6, experiment) algo.setMaximumOuterSampling(2500) algo.setBlockSize(4) algo.setMaximumCoefficientOfVariation(-1.0) algo.run() .. GENERATED FROM PYTHON SOURCE LINES 192-193 We retrieve the results and display the approximate probability and a confidence interval : .. GENERATED FROM PYTHON SOURCE LINES 193-200 .. code-block:: default 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)) .. rst-class:: sphx-glr-script-out Out: .. code-block:: none Probability of e6 through MC : 0.7474 Confidence interval MC : [0.7379, 0.7569] .. GENERATED FROM PYTHON SOURCE LINES 201-206 Usage with SystemFORM --------------------- The :class:`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). .. GENERATED FROM PYTHON SOURCE LINES 208-211 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. .. GENERATED FROM PYTHON SOURCE LINES 211-224 .. code-block:: default 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) .. GENERATED FROM PYTHON SOURCE LINES 225-226 As usual we create a :class:`~openturns.RandomVector` out of the input distribution. .. GENERATED FROM PYTHON SOURCE LINES 226-228 .. code-block:: default X = ot.RandomVector(dist) .. GENERATED FROM PYTHON SOURCE LINES 229-230 We define the leaf events thanks to :class:`~openturns.SymbolicFunction`. .. GENERATED FROM PYTHON SOURCE LINES 230-237 .. code-block:: default 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) .. GENERATED FROM PYTHON SOURCE LINES 238-239 We consider a system event in disjunctive normal form (union of intersections): .. GENERATED FROM PYTHON SOURCE LINES 239-241 .. code-block:: default event = ot.UnionEvent([ot.IntersectionEvent([e0, e3, e4]), ot.IntersectionEvent([e2, e3, e4])]) .. GENERATED FROM PYTHON SOURCE LINES 242-243 We can estimate the probability of the event with basic sampling. .. GENERATED FROM PYTHON SOURCE LINES 243-245 .. code-block:: default print("Probability of the event : %.4f"%event.getSample(10000).computeMean()[0]) .. rst-class:: sphx-glr-script-out Out: .. code-block:: none Probability of the event : 0.0769 .. GENERATED FROM PYTHON SOURCE LINES 246-247 We can also run a :class:`~openturns.systemFORM` algorithm to estimate the probability differently. .. GENERATED FROM PYTHON SOURCE LINES 249-250 We first set up a solver to find the design point. .. GENERATED FROM PYTHON SOURCE LINES 250-257 .. code-block:: default 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) .. GENERATED FROM PYTHON SOURCE LINES 258-259 We build the :class:`~openturns.SystemFORM` algorithm from the solver, the event and a starting point (here the mean) and then run the algorithm. .. GENERATED FROM PYTHON SOURCE LINES 259-262 .. code-block:: default algo = ot.SystemFORM(solver, event, mean) algo.run() .. GENERATED FROM PYTHON SOURCE LINES 263-264 We store the result and display the probability. .. GENERATED FROM PYTHON SOURCE LINES 264-268 .. code-block:: default result = algo.getResult() prbSystemFORM = result.getEventProbability() print("Probability of the event (SystemFORM) : %.4f"%prbSystemFORM) .. rst-class:: sphx-glr-script-out Out: .. code-block:: none Probability of the event (SystemFORM) : 0.0788 .. GENERATED FROM PYTHON SOURCE LINES 269-270 Display all figures .. GENERATED FROM PYTHON SOURCE LINES 270-271 .. code-block:: default plt.show() .. rst-class:: sphx-glr-timing **Total running time of the script:** ( 0 minutes 1.143 seconds) .. _sphx_glr_download_auto_reliability_sensitivity_reliability_plot_event_system.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_event_system.py ` .. container:: sphx-glr-download sphx-glr-download-jupyter :download:`Download Jupyter notebook: plot_event_system.ipynb ` .. only:: html .. rst-class:: sphx-glr-signature `Gallery generated by Sphinx-Gallery `_