.. 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 :ref:`Go to the end ` 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-14 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 14-21 .. code-block:: Python 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 22-36 **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 39-40 We consider a bivariate standard Gaussian random vector `X = (X_1, X_2)`. .. GENERATED FROM PYTHON SOURCE LINES 40-45 .. code-block:: Python dim = 2 distribution = ot.Normal(dim) X = ot.RandomVector(distribution) .. GENERATED FROM PYTHON SOURCE LINES 46-54 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 56-57 We consider three functions `f1`, `f2` and `f3` : .. GENERATED FROM PYTHON SOURCE LINES 57-61 .. code-block:: Python f1 = ot.SymbolicFunction(["x0", "x1"], ["x0"]) f2 = ot.SymbolicFunction(["x0", "x1"], ["x1"]) f3 = ot.SymbolicFunction(["x0", "x1"], ["x0+x1"]) .. GENERATED FROM PYTHON SOURCE LINES 62-63 We build :class:`~openturns.CompositeRandomVector` from these functions and the initial distribution. .. GENERATED FROM PYTHON SOURCE LINES 63-67 .. code-block:: Python Y1 = ot.CompositeRandomVector(f1, X) Y2 = ot.CompositeRandomVector(f2, X) Y3 = ot.CompositeRandomVector(f3, X) .. GENERATED FROM PYTHON SOURCE LINES 68-69 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 69-73 .. code-block:: Python 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 74-75 The restriction of the domain :math:`E_1` to :math:`[-4,4] \times [-4, 4]` is the grey area. .. GENERATED FROM PYTHON SOURCE LINES 75-87 .. code-block:: Python 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 88-89 The restriction of the domain :math:`E_2` to :math:`[-4,4] \times [-4, 4]` is the grey area. .. GENERATED FROM PYTHON SOURCE LINES 89-100 .. code-block:: Python 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 101-102 The restriction of the domain :math:`E_3` to :math:`[-4,4] \times [-4, 4]` is the grey area. .. GENERATED FROM PYTHON SOURCE LINES 102-113 .. code-block:: Python 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:: Python 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-136 .. code-block:: Python 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 137-138 The probability of that event is :math:`P_{E_4} = 1/4`. A basic estimator is: .. GENERATED FROM PYTHON SOURCE LINES 138-140 .. code-block:: Python print("Probability of e4 : %.4f" % e4.getSample(10000).computeMean()[0]) .. rst-class:: sphx-glr-script-out .. code-block:: none Probability of e4 : 0.2499 .. GENERATED FROM PYTHON SOURCE LINES 141-142 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 142-144 .. code-block:: Python e5 = ot.UnionEvent([e1, e2]) .. GENERATED FROM PYTHON SOURCE LINES 145-146 The restriction of the domain :math:`E_5` to :math:`[-4,4] \times [-4, 4]` is the grey area. .. GENERATED FROM PYTHON SOURCE LINES 146-163 .. code-block:: Python 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 164-165 The probability of that event is :math:`P_{E_5} = 3/4`. A basic estimator is: .. GENERATED FROM PYTHON SOURCE LINES 165-167 .. code-block:: Python print("Probability of e5 : %.4f" % e5.getSample(10000).computeMean()[0]) .. rst-class:: sphx-glr-script-out .. code-block:: none Probability of e5 : 0.7543 .. GENERATED FROM PYTHON SOURCE LINES 168-169 It supports recursion. Let's define :math:`E_6 = E_1 \bigcup (E_2 \bigcap E_3)`. .. GENERATED FROM PYTHON SOURCE LINES 169-172 .. code-block:: Python e6 = ot.UnionEvent([e1, ot.IntersectionEvent([e2, e3])]) .. GENERATED FROM PYTHON SOURCE LINES 173-174 First we draw the domain :math:`E_6 = E_1 \bigcup (E_2 \bigcap E_3)` : .. GENERATED FROM PYTHON SOURCE LINES 174-187 .. code-block:: Python 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 188-190 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 190-193 .. code-block:: Python print("Probability of e6 : %.4f" % e6.getSample(10000).computeMean()[0]) .. rst-class:: sphx-glr-script-out .. code-block:: none Probability of e6 : 0.7436 .. GENERATED FROM PYTHON SOURCE LINES 194-198 Usage with a Monte-Carlo algorithm ---------------------------------- Of course, we can use simulation algorithms with this kind of events. .. GENERATED FROM PYTHON SOURCE LINES 200-201 We set up a :class:`~openturns.MonteCarloExperiment` and a :class:`~openturns.ProbabilitySimulationAlgorithm` on the event :math:`E_6`. .. GENERATED FROM PYTHON SOURCE LINES 201-208 .. code-block:: Python 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 209-210 We retrieve the results and display the approximate probability and a confidence interval : .. GENERATED FROM PYTHON SOURCE LINES 210-217 .. code-block:: Python 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 .. code-block:: none Probability of e6 through MC : 0.7436 Confidence interval MC : [0.7340, 0.7532] .. GENERATED FROM PYTHON SOURCE LINES 218-223 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 225-228 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 228-241 .. code-block:: Python 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 242-243 As usual we create a :class:`~openturns.RandomVector` out of the input distribution. .. GENERATED FROM PYTHON SOURCE LINES 243-245 .. code-block:: Python X = ot.RandomVector(dist) .. GENERATED FROM PYTHON SOURCE LINES 246-247 We define the leaf events thanks to :class:`~openturns.SymbolicFunction`. .. GENERATED FROM PYTHON SOURCE LINES 247-274 .. code-block:: Python 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 275-276 We consider a system event in disjunctive normal form (union of intersections): .. GENERATED FROM PYTHON SOURCE LINES 276-280 .. code-block:: Python event = ot.UnionEvent( [ot.IntersectionEvent([e0, e3, e4]), ot.IntersectionEvent([e2, e3, e4])] ) .. GENERATED FROM PYTHON SOURCE LINES 281-282 We can estimate the probability of the event with basic sampling. .. GENERATED FROM PYTHON SOURCE LINES 282-284 .. code-block:: Python print("Probability of the event : %.4f" % event.getSample(10000).computeMean()[0]) .. rst-class:: sphx-glr-script-out .. code-block:: none Probability of the event : 0.0757 .. GENERATED FROM PYTHON SOURCE LINES 285-286 We can also run a :class:`~openturns.systemFORM` algorithm to estimate the probability differently. .. GENERATED FROM PYTHON SOURCE LINES 288-289 We first set up a solver to find the design point. .. GENERATED FROM PYTHON SOURCE LINES 289-296 .. code-block:: Python 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 297-298 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 298-301 .. code-block:: Python algo = ot.SystemFORM(solver, event, mean) algo.run() .. GENERATED FROM PYTHON SOURCE LINES 302-303 We store the result and display the probability. .. GENERATED FROM PYTHON SOURCE LINES 303-307 .. code-block:: Python result = algo.getResult() prbSystemFORM = result.getEventProbability() print("Probability of the event (SystemFORM) : %.4f" % prbSystemFORM) .. rst-class:: sphx-glr-script-out .. code-block:: none Probability of the event (SystemFORM) : 0.0788 .. GENERATED FROM PYTHON SOURCE LINES 308-309 Display all figures .. GENERATED FROM PYTHON SOURCE LINES 309-310 .. code-block:: Python plt.show() .. _sphx_glr_download_auto_reliability_sensitivity_reliability_plot_event_system.py: .. only:: html .. container:: sphx-glr-footer sphx-glr-footer-example .. container:: sphx-glr-download sphx-glr-download-jupyter :download:`Download Jupyter notebook: plot_event_system.ipynb ` .. container:: sphx-glr-download sphx-glr-download-python :download:`Download Python source code: plot_event_system.py `