.. 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 7-12 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 16-18 Definition of the model ----------------------- .. GENERATED FROM PYTHON SOURCE LINES 20-27 .. 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 28-29 We load the model from the usecases module : .. GENERATED FROM PYTHON SOURCE LINES 29-31 .. code-block:: Python sm = stressed_beam.AxialStressedBeam() .. GENERATED FROM PYTHON SOURCE LINES 32-33 The limit state function is defined as a symbolic function in the `model` parameter of the `AxialStressedBeam` data class: .. GENERATED FROM PYTHON SOURCE LINES 33-35 .. code-block:: Python limitStateFunction = sm.model .. GENERATED FROM PYTHON SOURCE LINES 36-37 Before using the function within an algorithm, we check that the limit state function is correctly evaluated. .. GENERATED FROM PYTHON SOURCE LINES 39-43 .. 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 44-46 Probabilistic model ------------------- .. GENERATED FROM PYTHON SOURCE LINES 48-49 We load the first marginal, a univariate `LogNormal` distribution, parameterized by its mean and standard deviation: .. GENERATED FROM PYTHON SOURCE LINES 49-51 .. code-block:: Python R = sm.distribution_R .. GENERATED FROM PYTHON SOURCE LINES 52-53 We draw the PDF of the first marginal. .. GENERATED FROM PYTHON SOURCE LINES 53-56 .. 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 57-58 Our second marginal is a `Normal` univariate distribution: .. GENERATED FROM PYTHON SOURCE LINES 58-60 .. code-block:: Python F = sm.distribution_F .. GENERATED FROM PYTHON SOURCE LINES 61-62 We draw the PDF of the second marginal. .. GENERATED FROM PYTHON SOURCE LINES 62-65 .. 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 66-70 In order to create the input distribution, we use the `JointDistribution` 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 72-74 .. code-block:: Python myDistribution = sm.distribution .. GENERATED FROM PYTHON SOURCE LINES 75-76 We create a `RandomVector` from the `Distribution`, which will then be fed to the limit state function. .. GENERATED FROM PYTHON SOURCE LINES 78-80 .. code-block:: Python inputRandomVector = ot.RandomVector(myDistribution) .. GENERATED FROM PYTHON SOURCE LINES 81-82 Finally we create a `CompositeRandomVector` by associating the limit state function with the input random vector. .. GENERATED FROM PYTHON SOURCE LINES 84-86 .. code-block:: Python outputRandomVector = ot.CompositeRandomVector(limitStateFunction, inputRandomVector) .. GENERATED FROM PYTHON SOURCE LINES 87-89 Exact computation ----------------- .. GENERATED FROM PYTHON SOURCE LINES 91-92 The simplest method is to perform an exact computation based on the arithmetic of distributions. .. GENERATED FROM PYTHON SOURCE LINES 94-96 .. code-block:: Python D = 0.02 .. GENERATED FROM PYTHON SOURCE LINES 97-99 .. code-block:: Python G = R - F / (D**2 / 4 * np.pi) .. GENERATED FROM PYTHON SOURCE LINES 100-102 .. code-block:: Python G.computeCDF(0.0) .. rst-class:: sphx-glr-script-out .. code-block:: none 0.029198194624830504 .. GENERATED FROM PYTHON SOURCE LINES 103-104 This is the exact result from the description of this example. .. GENERATED FROM PYTHON SOURCE LINES 106-108 Distribution of the output -------------------------- .. GENERATED FROM PYTHON SOURCE LINES 110-111 Plot the distribution of the output. .. GENERATED FROM PYTHON SOURCE LINES 113-118 .. 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 119-121 Estimate the probability with Monte-Carlo ----------------------------------------- .. GENERATED FROM PYTHON SOURCE LINES 123-124 We first create a `ThresholdEvent` based on the output `RandomVector`, the operator and the threshold. .. GENERATED FROM PYTHON SOURCE LINES 126-128 .. code-block:: Python myEvent = ot.ThresholdEvent(outputRandomVector, ot.Less(), 0.0) .. GENERATED FROM PYTHON SOURCE LINES 129-131 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 133-142 .. 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 143-145 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 147-149 .. code-block:: Python initialNumberOfCall = limitStateFunction.getEvaluationCallsNumber() .. GENERATED FROM PYTHON SOURCE LINES 150-152 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 154-156 .. code-block:: Python algoMC.run() .. GENERATED FROM PYTHON SOURCE LINES 157-158 We can then get the results of the algorithm. .. GENERATED FROM PYTHON SOURCE LINES 160-169 .. 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 = 14897 Pf = 0.026179767738470774 CV = 0.04996974038872697 .. GENERATED FROM PYTHON SOURCE LINES 170-173 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 175-179 .. 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 180-184 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 186-188 .. code-block:: Python alpha = 0.05 .. GENERATED FROM PYTHON SOURCE LINES 189-195 .. 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.023616,0.028744] .. GENERATED FROM PYTHON SOURCE LINES 196-197 This interval is consistent with the exact probability :math:`P_f=0.02920`. .. GENERATED FROM PYTHON SOURCE LINES 199-237 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 ` .. container:: sphx-glr-download sphx-glr-download-zip :download:`Download zipped: plot_axial_stressed_beam_quickstart.zip `