.. 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-17 .. code-block:: default 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 18-32 **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 35-36 We consider a bivariate standard Gaussian random vector `X = (X_1, X_2)`. .. GENERATED FROM PYTHON SOURCE LINES 36-41 .. code-block:: default dim = 2 distribution = ot.Normal(dim) X = ot.RandomVector(distribution) .. GENERATED FROM PYTHON SOURCE LINES 42-50 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 52-53 We consider three functions `f1`, `f2` and `f3` : .. GENERATED FROM PYTHON SOURCE LINES 53-57 .. 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 58-59 We build :class:`~openturns.CompositeRandomVector` from these functions and the initial distribution. .. GENERATED FROM PYTHON SOURCE LINES 59-63 .. code-block:: default Y1 = ot.CompositeRandomVector(f1, X) Y2 = ot.CompositeRandomVector(f2, X) Y3 = ot.CompositeRandomVector(f3, X) .. GENERATED FROM PYTHON SOURCE LINES 64-65 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 65-69 .. 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 70-71 The restriction of the domain :math:`E_1` to :math:`[-4,4] \times [-4, 4]` is the grey area. .. GENERATED FROM PYTHON SOURCE LINES 71-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-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 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-98 .. 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 99-100 The restriction of the domain :math:`E_3` to :math:`[-4,4] \times [-4, 4]` is the grey area. .. GENERATED FROM PYTHON SOURCE LINES 100-112 .. 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 113-114 We can define the intersection :math:`E_4 = E_1 \bigcap E_2`: that is the upper left quadrant. .. GENERATED FROM PYTHON SOURCE LINES 114-116 .. code-block:: default e4 = ot.IntersectionEvent([e1, e2]) .. GENERATED FROM PYTHON SOURCE LINES 117-118 The restriction of the domain :math:`E_4` to :math:`[-4,4] \times [-4, 4]` is the grey area. .. GENERATED FROM PYTHON SOURCE LINES 118-130 .. 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 131-132 The probability of that event is :math:`P_{E_4} = 1/4`. A basic estimator is: .. GENERATED FROM PYTHON SOURCE LINES 132-134 .. 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.2413 .. GENERATED FROM PYTHON SOURCE LINES 135-136 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 136-138 .. code-block:: default e5 = ot.UnionEvent([e1, e2]) .. GENERATED FROM PYTHON SOURCE LINES 139-140 The restriction of the domain :math:`E_5` to :math:`[-4,4] \times [-4, 4]` is the grey area. .. GENERATED FROM PYTHON SOURCE LINES 140-152 .. 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 153-154 The probability of that event is :math:`P_{E_5} = 3/4`. A basic estimator is: .. GENERATED FROM PYTHON SOURCE LINES 154-156 .. 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.7494 .. GENERATED FROM PYTHON SOURCE LINES 157-158 It supports recursion. Let's define :math:`E_6 = E_1 \bigcup (E_2 \bigcap E_3)`. .. GENERATED FROM PYTHON SOURCE LINES 158-161 .. code-block:: default e6 = ot.UnionEvent([e1, ot.IntersectionEvent([e2, e3])]) .. GENERATED FROM PYTHON SOURCE LINES 162-163 First we draw the domain :math:`E_6 = E_1 \bigcup (E_2 \bigcap E_3)` : .. GENERATED FROM PYTHON SOURCE LINES 163-175 .. 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 176-178 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 178-181 .. 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.7511 .. GENERATED FROM PYTHON SOURCE LINES 182-186 Usage with a Monte-Carlo algorithm ---------------------------------- Of course, we can use simulation algorithms with this kind of events. .. GENERATED FROM PYTHON SOURCE LINES 188-189 We set up a :class:`~openturns.MonteCarloExperiment` and a :class:`~openturns.ProbabilitySimulationAlgorithm` on the event :math:`E_6`. .. GENERATED FROM PYTHON SOURCE LINES 189-196 .. 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 197-198 We retrieve the results and display the approximate probability and a confidence interval : .. GENERATED FROM PYTHON SOURCE LINES 198-205 .. 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.7508 Confidence interval MC : [0.7413, 0.7603] .. GENERATED FROM PYTHON SOURCE LINES 206-211 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 213-216 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 216-229 .. 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 230-231 As usual we create a :class:`~openturns.RandomVector` out of the input distribution. .. GENERATED FROM PYTHON SOURCE LINES 231-233 .. code-block:: default X = ot.RandomVector(dist) .. GENERATED FROM PYTHON SOURCE LINES 234-235 We define the leaf events thanks to :class:`~openturns.SymbolicFunction`. .. GENERATED FROM PYTHON SOURCE LINES 235-247 .. 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 248-249 We consider a system event in disjunctive normal form (union of intersections): .. GENERATED FROM PYTHON SOURCE LINES 249-252 .. code-block:: default event = ot.UnionEvent([ot.IntersectionEvent([e0, e3, e4]), ot.IntersectionEvent([e2, e3, e4])]) .. GENERATED FROM PYTHON SOURCE LINES 253-254 We can estimate the probability of the event with basic sampling. .. GENERATED FROM PYTHON SOURCE LINES 254-257 .. 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.0789 .. GENERATED FROM PYTHON SOURCE LINES 258-259 We can also run a :class:`~openturns.systemFORM` algorithm to estimate the probability differently. .. GENERATED FROM PYTHON SOURCE LINES 261-262 We first set up a solver to find the design point. .. GENERATED FROM PYTHON SOURCE LINES 262-269 .. 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 270-271 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 271-274 .. code-block:: default algo = ot.SystemFORM(solver, event, mean) algo.run() .. GENERATED FROM PYTHON SOURCE LINES 275-276 We store the result and display the probability. .. GENERATED FROM PYTHON SOURCE LINES 276-280 .. 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 281-282 Display all figures .. GENERATED FROM PYTHON SOURCE LINES 282-283 .. code-block:: default plt.show() .. rst-class:: sphx-glr-timing **Total running time of the script:** ( 0 minutes 0.722 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 `_