.. 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 Click :ref:`here ` 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-8 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 12-14 Definition of the model ----------------------- .. GENERATED FROM PYTHON SOURCE LINES 16-23 .. code-block:: default from openturns.usecases import stressed_beam as stressed_beam import openturns as ot import numpy as np import openturns.viewer as viewer from matplotlib import pylab as plt ot.Log.Show(ot.Log.NONE) .. GENERATED FROM PYTHON SOURCE LINES 24-25 We load the model from the usecases module : .. GENERATED FROM PYTHON SOURCE LINES 25-27 .. code-block:: default sm = stressed_beam.AxialStressedBeam() .. GENERATED FROM PYTHON SOURCE LINES 28-29 The limit state function is defined as a symbolic function in the `model` parameter of the `AxialStressedBeam` data class : .. GENERATED FROM PYTHON SOURCE LINES 29-31 .. code-block:: default limitStateFunction = sm.model .. GENERATED FROM PYTHON SOURCE LINES 32-33 Before using the function within an algorithm, we check that the limit state function is correctly evaluated. .. GENERATED FROM PYTHON SOURCE LINES 35-39 .. code-block:: default x = [3.e6, 750.] print('x=', x) print('G(x)=', limitStateFunction(x)) .. rst-class:: sphx-glr-script-out Out: .. code-block:: none x= [3000000.0, 750.0] G(x)= [612676] .. GENERATED FROM PYTHON SOURCE LINES 40-42 Probabilistic model ------------------- .. GENERATED FROM PYTHON SOURCE LINES 44-45 We load the first marginal, a univariate `LogNormal` distribution, parameterized by its mean and standard deviation : .. GENERATED FROM PYTHON SOURCE LINES 45-47 .. code-block:: default R = sm.distribution_R .. GENERATED FROM PYTHON SOURCE LINES 48-49 We draw the PDF of the first marginal. .. GENERATED FROM PYTHON SOURCE LINES 49-52 .. code-block:: default 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 53-54 Our second marginal is a `Normal` univariate distribution : .. GENERATED FROM PYTHON SOURCE LINES 54-56 .. code-block:: default F = sm.distribution_F .. GENERATED FROM PYTHON SOURCE LINES 57-58 We draw the PDF of the second marginal. .. GENERATED FROM PYTHON SOURCE LINES 58-61 .. code-block:: default 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 62-63 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 65-67 .. code-block:: default myDistribution = sm.distribution .. GENERATED FROM PYTHON SOURCE LINES 68-69 We create a `RandomVector` from the `Distribution`, which will then be fed to the limit state function. .. GENERATED FROM PYTHON SOURCE LINES 71-73 .. code-block:: default inputRandomVector = ot.RandomVector(myDistribution) .. GENERATED FROM PYTHON SOURCE LINES 74-75 Finally we create a `CompositeRandomVector` by associating the limit state function with the input random vector. .. GENERATED FROM PYTHON SOURCE LINES 77-80 .. code-block:: default outputRandomVector = ot.CompositeRandomVector( limitStateFunction, inputRandomVector) .. GENERATED FROM PYTHON SOURCE LINES 81-83 Exact computation ----------------- .. GENERATED FROM PYTHON SOURCE LINES 85-86 The simplest method is to perform an exact computation based on the arithmetic of distributions. .. GENERATED FROM PYTHON SOURCE LINES 88-90 .. code-block:: default D = 0.02 .. GENERATED FROM PYTHON SOURCE LINES 91-93 .. code-block:: default G = R-F/(D**2/4 * np.pi) .. GENERATED FROM PYTHON SOURCE LINES 94-96 .. code-block:: default G.computeCDF(0.) .. rst-class:: sphx-glr-script-out Out: .. code-block:: none 0.029198194624830504 .. GENERATED FROM PYTHON SOURCE LINES 97-98 This is the exact result from the description of this example. .. GENERATED FROM PYTHON SOURCE LINES 100-102 Distribution of the output -------------------------- .. GENERATED FROM PYTHON SOURCE LINES 104-105 Plot the distribution of the output. .. GENERATED FROM PYTHON SOURCE LINES 107-112 .. code-block:: default 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 113-115 Estimate the probability with Monte-Carlo ----------------------------------------- .. GENERATED FROM PYTHON SOURCE LINES 117-118 We first create a `ThresholdEvent` based on the output `RandomVector`, the operator and the threshold. .. GENERATED FROM PYTHON SOURCE LINES 120-122 .. code-block:: default myEvent = ot.ThresholdEvent(outputRandomVector, ot.Less(), 0.0) .. GENERATED FROM PYTHON SOURCE LINES 123-124 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 126-135 .. code-block:: default 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 136-137 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 139-141 .. code-block:: default initialNumberOfCall = limitStateFunction.getEvaluationCallsNumber() .. GENERATED FROM PYTHON SOURCE LINES 142-143 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 145-147 .. code-block:: default algoMC.run() .. GENERATED FROM PYTHON SOURCE LINES 148-149 We can then get the results of the algorithm. .. GENERATED FROM PYTHON SOURCE LINES 151-159 .. code-block:: default 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 Out: .. code-block:: none Number of calls to the limit state = 11897 Pf = 0.032529209044296944 CV = 0.049999245238636636 .. GENERATED FROM PYTHON SOURCE LINES 160-161 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 163-167 .. code-block:: default 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 168-171 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 173-175 .. code-block:: default alpha = 0.05 .. GENERATED FROM PYTHON SOURCE LINES 176-180 .. code-block:: default 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 Out: .. code-block:: none 95.00% confidence interval = [0.029341,0.035717] .. GENERATED FROM PYTHON SOURCE LINES 181-182 This interval is consistent with the exact probability :math:`P_f=0.02920`. .. GENERATED FROM PYTHON SOURCE LINES 184-221 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`. .. rst-class:: sphx-glr-timing **Total running time of the script:** ( 0 minutes 0.650 seconds) .. _sphx_glr_download_auto_reliability_sensitivity_reliability_plot_axial_stressed_beam_quickstart.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_axial_stressed_beam_quickstart.py ` .. container:: sphx-glr-download sphx-glr-download-jupyter :download:`Download Jupyter notebook: plot_axial_stressed_beam_quickstart.ipynb ` .. only:: html .. rst-class:: sphx-glr-signature `Gallery generated by Sphinx-Gallery `_