.. 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-85 .. 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-sg:: /auto_reliability_sensitivity/reliability/images/sphx_glr_plot_event_system_001.png :alt: Representation of the event $E_1$ :srcset: /auto_reliability_sensitivity/reliability/images/sphx_glr_plot_event_system_001.png :class: sphx-glr-single-img .. GENERATED FROM PYTHON SOURCE LINES 86-87 The restriction of the domain :math:`E_2` to :math:`[-4,4] \times [-4, 4]` is the grey area. .. GENERATED FROM PYTHON SOURCE LINES 87-99 .. 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-sg:: /auto_reliability_sensitivity/reliability/images/sphx_glr_plot_event_system_002.png :alt: Representation of the event $E_2$ :srcset: /auto_reliability_sensitivity/reliability/images/sphx_glr_plot_event_system_002.png :class: sphx-glr-single-img .. GENERATED FROM PYTHON SOURCE LINES 100-101 The restriction of the domain :math:`E_3` to :math:`[-4,4] \times [-4, 4]` is the grey area. .. GENERATED FROM PYTHON SOURCE LINES 101-113 .. 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-sg:: /auto_reliability_sensitivity/reliability/images/sphx_glr_plot_event_system_003.png :alt: Representation of the event $E_3$ :srcset: /auto_reliability_sensitivity/reliability/images/sphx_glr_plot_event_system_003.png :class: sphx-glr-single-img .. GENERATED FROM PYTHON SOURCE LINES 114-115 We can define the intersection :math:`E_4 = E_1 \bigcap E_2`: that is the upper left quadrant. .. GENERATED FROM PYTHON SOURCE LINES 115-117 .. code-block:: default e4 = ot.IntersectionEvent([e1, e2]) .. GENERATED FROM PYTHON SOURCE LINES 118-119 The restriction of the domain :math:`E_4` to :math:`[-4,4] \times [-4, 4]` is the grey area. .. GENERATED FROM PYTHON SOURCE LINES 119-131 .. 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-sg:: /auto_reliability_sensitivity/reliability/images/sphx_glr_plot_event_system_004.png :alt: Representation of the event $E_4 = E_1 \bigcap E_2$ :srcset: /auto_reliability_sensitivity/reliability/images/sphx_glr_plot_event_system_004.png :class: sphx-glr-single-img .. GENERATED FROM PYTHON SOURCE LINES 132-133 The probability of that event is :math:`P_{E_4} = 1/4`. A basic estimator is: .. GENERATED FROM PYTHON SOURCE LINES 133-135 .. 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.2504 .. GENERATED FROM PYTHON SOURCE LINES 136-137 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 137-139 .. code-block:: default e5 = ot.UnionEvent([e1, e2]) .. GENERATED FROM PYTHON SOURCE LINES 140-141 The restriction of the domain :math:`E_5` to :math:`[-4,4] \times [-4, 4]` is the grey area. .. GENERATED FROM PYTHON SOURCE LINES 141-153 .. 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-sg:: /auto_reliability_sensitivity/reliability/images/sphx_glr_plot_event_system_005.png :alt: Representation of the event $E_5 = E_1 \bigcup E_2$ :srcset: /auto_reliability_sensitivity/reliability/images/sphx_glr_plot_event_system_005.png :class: sphx-glr-single-img .. GENERATED FROM PYTHON SOURCE LINES 154-155 The probability of that event is :math:`P_{E_5} = 3/4`. A basic estimator is: .. GENERATED FROM PYTHON SOURCE LINES 155-157 .. 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.7489 .. GENERATED FROM PYTHON SOURCE LINES 158-159 It supports recursion. Let's define :math:`E_6 = E_1 \bigcup (E_2 \bigcap E_3)`. .. GENERATED FROM PYTHON SOURCE LINES 159-162 .. code-block:: default e6 = ot.UnionEvent([e1, ot.IntersectionEvent([e2, e3])]) .. GENERATED FROM PYTHON SOURCE LINES 163-164 First we draw the domain :math:`E_6 = E_1 \bigcup (E_2 \bigcap E_3)` : .. GENERATED FROM PYTHON SOURCE LINES 164-176 .. 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-sg:: /auto_reliability_sensitivity/reliability/images/sphx_glr_plot_event_system_006.png :alt: Representation of the event $E_2 \bigcap E_3 $ :srcset: /auto_reliability_sensitivity/reliability/images/sphx_glr_plot_event_system_006.png :class: sphx-glr-single-img .. GENERATED FROM PYTHON SOURCE LINES 177-179 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 179-182 .. 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.7510 .. GENERATED FROM PYTHON SOURCE LINES 183-187 Usage with a Monte-Carlo algorithm ---------------------------------- Of course, we can use simulation algorithms with this kind of events. .. GENERATED FROM PYTHON SOURCE LINES 189-190 We set up a :class:`~openturns.MonteCarloExperiment` and a :class:`~openturns.ProbabilitySimulationAlgorithm` on the event :math:`E_6`. .. GENERATED FROM PYTHON SOURCE LINES 190-197 .. 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 198-199 We retrieve the results and display the approximate probability and a confidence interval : .. GENERATED FROM PYTHON SOURCE LINES 199-206 .. 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.7483 Confidence interval MC : [0.7388, 0.7578] .. GENERATED FROM PYTHON SOURCE LINES 207-212 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 214-217 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 217-230 .. 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 231-232 As usual we create a :class:`~openturns.RandomVector` out of the input distribution. .. GENERATED FROM PYTHON SOURCE LINES 232-234 .. code-block:: default X = ot.RandomVector(dist) .. GENERATED FROM PYTHON SOURCE LINES 235-236 We define the leaf events thanks to :class:`~openturns.SymbolicFunction`. .. GENERATED FROM PYTHON SOURCE LINES 236-248 .. 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 249-250 We consider a system event in disjunctive normal form (union of intersections): .. GENERATED FROM PYTHON SOURCE LINES 250-253 .. code-block:: default event = ot.UnionEvent([ot.IntersectionEvent([e0, e3, e4]), ot.IntersectionEvent([e2, e3, e4])]) .. GENERATED FROM PYTHON SOURCE LINES 254-255 We can estimate the probability of the event with basic sampling. .. GENERATED FROM PYTHON SOURCE LINES 255-258 .. 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.0788 .. GENERATED FROM PYTHON SOURCE LINES 259-260 We can also run a :class:`~openturns.systemFORM` algorithm to estimate the probability differently. .. GENERATED FROM PYTHON SOURCE LINES 262-263 We first set up a solver to find the design point. .. GENERATED FROM PYTHON SOURCE LINES 263-270 .. 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 271-272 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 272-275 .. code-block:: default algo = ot.SystemFORM(solver, event, mean) algo.run() .. GENERATED FROM PYTHON SOURCE LINES 276-277 We store the result and display the probability. .. GENERATED FROM PYTHON SOURCE LINES 277-281 .. 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 282-283 Display all figures .. GENERATED FROM PYTHON SOURCE LINES 283-284 .. code-block:: default plt.show() .. rst-class:: sphx-glr-timing **Total running time of the script:** ( 0 minutes 1.284 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 `_