.. DO NOT EDIT. .. THIS FILE WAS AUTOMATICALLY GENERATED BY SPHINX-GALLERY. .. TO MAKE CHANGES, EDIT THE SOURCE PYTHON FILE: .. "auto_reliability_sensitivity/reliability/plot_axial_stressed_beam_quickstart.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_axial_stressed_beam_quickstart.py: Estimate a probability with Monte-Carlo on axial stressed beam: a quick start guide to reliability ================================================================================================== .. GENERATED FROM PYTHON SOURCE LINES 6-11 The goal of this example is to show a simple practical example of probability estimation in a reliability study with the `ProbabilitySimulationAlgorithm` class. The `ThresholdEvent` is used to define the event. We use the Monte-Carlo method thanks to the `MonteCarloExperiment` class to estimate this probability and its confidence interval. We use the :ref:`axial stressed beam ` model as an illustrative example. .. GENERATED FROM PYTHON SOURCE LINES 15-17 Definition of the model ----------------------- .. GENERATED FROM PYTHON SOURCE LINES 19-26 .. code-block:: Python from openturns.usecases import stressed_beam import openturns as ot import numpy as np import openturns.viewer as viewer ot.Log.Show(ot.Log.NONE) .. GENERATED FROM PYTHON SOURCE LINES 27-28 We load the model from the usecases module : .. GENERATED FROM PYTHON SOURCE LINES 28-30 .. code-block:: Python sm = stressed_beam.AxialStressedBeam() .. GENERATED FROM PYTHON SOURCE LINES 31-32 The limit state function is defined as a symbolic function in the `model` parameter of the `AxialStressedBeam` data class: .. GENERATED FROM PYTHON SOURCE LINES 32-34 .. code-block:: Python limitStateFunction = sm.model .. GENERATED FROM PYTHON SOURCE LINES 35-36 Before using the function within an algorithm, we check that the limit state function is correctly evaluated. .. GENERATED FROM PYTHON SOURCE LINES 38-42 .. code-block:: Python x = [3.0e6, 750.0] print("x=", x) print("G(x)=", limitStateFunction(x)) .. rst-class:: sphx-glr-script-out .. code-block:: none x= [3000000.0, 750.0] G(x)= [612676] .. GENERATED FROM PYTHON SOURCE LINES 43-45 Probabilistic model ------------------- .. GENERATED FROM PYTHON SOURCE LINES 47-48 We load the first marginal, a univariate `LogNormal` distribution, parameterized by its mean and standard deviation: .. GENERATED FROM PYTHON SOURCE LINES 48-50 .. code-block:: Python R = sm.distribution_R .. GENERATED FROM PYTHON SOURCE LINES 51-52 We draw the PDF of the first marginal. .. GENERATED FROM PYTHON SOURCE LINES 52-55 .. code-block:: Python graph = R.drawPDF() view = viewer.View(graph) .. image-sg:: /auto_reliability_sensitivity/reliability/images/sphx_glr_plot_axial_stressed_beam_quickstart_001.png :alt: plot axial stressed beam quickstart :srcset: /auto_reliability_sensitivity/reliability/images/sphx_glr_plot_axial_stressed_beam_quickstart_001.png :class: sphx-glr-single-img .. GENERATED FROM PYTHON SOURCE LINES 56-57 Our second marginal is a `Normal` univariate distribution: .. GENERATED FROM PYTHON SOURCE LINES 57-59 .. code-block:: Python F = sm.distribution_F .. GENERATED FROM PYTHON SOURCE LINES 60-61 We draw the PDF of the second marginal. .. GENERATED FROM PYTHON SOURCE LINES 61-64 .. code-block:: Python graph = F.drawPDF() view = viewer.View(graph) .. image-sg:: /auto_reliability_sensitivity/reliability/images/sphx_glr_plot_axial_stressed_beam_quickstart_002.png :alt: plot axial stressed beam quickstart :srcset: /auto_reliability_sensitivity/reliability/images/sphx_glr_plot_axial_stressed_beam_quickstart_002.png :class: sphx-glr-single-img .. GENERATED FROM PYTHON SOURCE LINES 65-69 In order to create the input distribution, we use the `ComposedDistribution` class which associates the distribution marginals and a copula. If no copula is supplied to the constructor, it selects the independent copula as default. That is implemented in the data class: .. GENERATED FROM PYTHON SOURCE LINES 71-73 .. code-block:: Python myDistribution = sm.distribution .. GENERATED FROM PYTHON SOURCE LINES 74-75 We create a `RandomVector` from the `Distribution`, which will then be fed to the limit state function. .. GENERATED FROM PYTHON SOURCE LINES 77-79 .. code-block:: Python inputRandomVector = ot.RandomVector(myDistribution) .. GENERATED FROM PYTHON SOURCE LINES 80-81 Finally we create a `CompositeRandomVector` by associating the limit state function with the input random vector. .. GENERATED FROM PYTHON SOURCE LINES 83-85 .. code-block:: Python outputRandomVector = ot.CompositeRandomVector(limitStateFunction, inputRandomVector) .. GENERATED FROM PYTHON SOURCE LINES 86-88 Exact computation ----------------- .. GENERATED FROM PYTHON SOURCE LINES 90-91 The simplest method is to perform an exact computation based on the arithmetic of distributions. .. GENERATED FROM PYTHON SOURCE LINES 93-95 .. code-block:: Python D = 0.02 .. GENERATED FROM PYTHON SOURCE LINES 96-98 .. code-block:: Python G = R - F / (D**2 / 4 * np.pi) .. GENERATED FROM PYTHON SOURCE LINES 99-101 .. code-block:: Python G.computeCDF(0.0) .. rst-class:: sphx-glr-script-out .. code-block:: none 0.029198194624830504 .. GENERATED FROM PYTHON SOURCE LINES 102-103 This is the exact result from the description of this example. .. GENERATED FROM PYTHON SOURCE LINES 105-107 Distribution of the output -------------------------- .. GENERATED FROM PYTHON SOURCE LINES 109-110 Plot the distribution of the output. .. GENERATED FROM PYTHON SOURCE LINES 112-117 .. code-block:: Python sampleSize = 500 sampleG = outputRandomVector.getSample(sampleSize) graph = ot.HistogramFactory().build(sampleG).drawPDF() view = viewer.View(graph) .. image-sg:: /auto_reliability_sensitivity/reliability/images/sphx_glr_plot_axial_stressed_beam_quickstart_003.png :alt: y0 PDF :srcset: /auto_reliability_sensitivity/reliability/images/sphx_glr_plot_axial_stressed_beam_quickstart_003.png :class: sphx-glr-single-img .. GENERATED FROM PYTHON SOURCE LINES 118-120 Estimate the probability with Monte-Carlo ----------------------------------------- .. GENERATED FROM PYTHON SOURCE LINES 122-123 We first create a `ThresholdEvent` based on the output `RandomVector`, the operator and the threshold. .. GENERATED FROM PYTHON SOURCE LINES 125-127 .. code-block:: Python myEvent = ot.ThresholdEvent(outputRandomVector, ot.Less(), 0.0) .. GENERATED FROM PYTHON SOURCE LINES 128-130 The `ProbabilitySimulationAlgorithm` is the main tool to estimate a probability. It is based on a specific design of experiments: in this example, we use the simplest of all, the `MonteCarloExperiment`. .. GENERATED FROM PYTHON SOURCE LINES 132-141 .. code-block:: Python maximumCoV = 0.05 # Coefficient of variation maximumNumberOfBlocks = 100000 experiment = ot.MonteCarloExperiment() algoMC = ot.ProbabilitySimulationAlgorithm(myEvent, experiment) algoMC.setMaximumOuterSampling(maximumNumberOfBlocks) algoMC.setBlockSize(1) algoMC.setMaximumCoefficientOfVariation(maximumCoV) .. GENERATED FROM PYTHON SOURCE LINES 142-144 In order to gather statistics about the algorithm, we get the initial number of function calls (this is not mandatory, but will prove to be convenient later in the study). .. GENERATED FROM PYTHON SOURCE LINES 146-148 .. code-block:: Python initialNumberOfCall = limitStateFunction.getEvaluationCallsNumber() .. GENERATED FROM PYTHON SOURCE LINES 149-151 Now comes the costly part: the `run` method performs the required simulations. The algorithm stops when the coefficient of variation of the probability estimate becomes lower than 0.5. .. GENERATED FROM PYTHON SOURCE LINES 153-155 .. code-block:: Python algoMC.run() .. GENERATED FROM PYTHON SOURCE LINES 156-157 We can then get the results of the algorithm. .. GENERATED FROM PYTHON SOURCE LINES 159-168 .. code-block:: Python result = algoMC.getResult() probability = result.getProbabilityEstimate() numberOfFunctionEvaluations = ( limitStateFunction.getEvaluationCallsNumber() - initialNumberOfCall ) print("Number of calls to the limit state =", numberOfFunctionEvaluations) print("Pf = ", probability) print("CV =", result.getCoefficientOfVariation()) .. rst-class:: sphx-glr-script-out .. code-block:: none Number of calls to the limit state = 13424 Pf = 0.028977949940405243 CV = 0.049961991982172875 .. GENERATED FROM PYTHON SOURCE LINES 169-172 The `drawProbabilityConvergence` method plots the probability estimate depending on the number of function evaluations. The order of convergence is :math:`O \left( 1/N^2 \right)` with :math:`N` being the number of function evaluations. This is why we use a logarithmic scale for the X axis of the graphics. .. GENERATED FROM PYTHON SOURCE LINES 174-178 .. code-block:: Python graph = algoMC.drawProbabilityConvergence() graph.setLogScale(ot.GraphImplementation.LOGX) view = viewer.View(graph) .. image-sg:: /auto_reliability_sensitivity/reliability/images/sphx_glr_plot_axial_stressed_beam_quickstart_004.png :alt: ProbabilitySimulationAlgorithm convergence graph at level 0.95 :srcset: /auto_reliability_sensitivity/reliability/images/sphx_glr_plot_axial_stressed_beam_quickstart_004.png :class: sphx-glr-single-img .. GENERATED FROM PYTHON SOURCE LINES 179-183 We see that the 95% confidence interval becomes smaller and smaller and stabilizes at the end of the simulation. In order to compute the confidence interval, we use the `getConfidenceLength` method, which returns the length of the interval. In order to compute the bounds of the interval, we divide this length by 2. .. GENERATED FROM PYTHON SOURCE LINES 185-187 .. code-block:: Python alpha = 0.05 .. GENERATED FROM PYTHON SOURCE LINES 188-194 .. code-block:: Python pflen = result.getConfidenceLength(1 - alpha) print( "%.2f%% confidence interval = [%f,%f]" % ((1 - alpha) * 100, probability - pflen / 2, probability + pflen / 2) ) .. rst-class:: sphx-glr-script-out .. code-block:: none 95.00% confidence interval = [0.026140,0.031816] .. GENERATED FROM PYTHON SOURCE LINES 195-196 This interval is consistent with the exact probability :math:`P_f=0.02920`. .. GENERATED FROM PYTHON SOURCE LINES 198-236 Appendix: derivation of the failure probability ----------------------------------------------- The failure probability is: .. math:: P_f = \text{Prob}(R-S \leq 0) = \int_{r-s \leq 0} f_{R, S}(r, s)drds where :math:`f_{R, S}` is the probability distribution function of the random vector :math:`(R,S)`. If R and S are independent, then: .. math:: f_{R, S}(r, s) = f_R(r) f_S(s) for any :math:`r,s\in\mathbb{R}`, where :math:`f_S` is the probability distribution function of the random variable :math:`S` and :math:`f_R` is the probability distribution function of the random variable :math:`R`. Therefore, .. math:: P_f = \int_{r-s \leq 0} f_R(r) f_S(s) drds. This implies: .. math:: P_f = \int_{-\infty}^{+\infty} \left(\int_{r \leq s} f_R(r) dr \right) f_S(s) ds. Therefore, .. math:: P_f = \int_{-\infty}^{+\infty}f_S(s)F_R(s)ds where :math:`F_R` is the cumulative distribution function of the random variable :math:`R`. .. _sphx_glr_download_auto_reliability_sensitivity_reliability_plot_axial_stressed_beam_quickstart.py: .. only:: html .. container:: sphx-glr-footer sphx-glr-footer-example .. container:: sphx-glr-download sphx-glr-download-jupyter :download:`Download Jupyter notebook: plot_axial_stressed_beam_quickstart.ipynb ` .. container:: sphx-glr-download sphx-glr-download-python :download:`Download Python source code: plot_axial_stressed_beam_quickstart.py `