.. only:: html

    .. note::
        :class: sphx-glr-download-link-note

        Click :ref:`here <sphx_glr_download_auto_reliability_sensitivity_reliability_plot_axial_stressed_beam_quickstart.py>`     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
==================================================================================================

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 <use-case-stressed-beam>` model as an illustrative example.

Definition of the model
-----------------------


.. code-block:: default

    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)








We load the model from the usecases module :


.. code-block:: default

    from openturns.usecases import stressed_beam as stressed_beam
    sm = stressed_beam.AxialStressedBeam()








The limit state function is defined as a symbolic function in the `model` parameter of the `AxialStressedBeam` data class : 


.. code-block:: default

    limitStateFunction = sm.model








Before using the function within an algorithm, we check that the limit state function is correctly evaluated.


.. 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]




Probabilistic model
-------------------

We load the first marginal, a univariate `LogNormal` distribution, parameterized by its mean and standard deviation :


.. code-block:: default

    R = sm.distribution_R








We draw the PDF of the first marginal.


.. code-block:: default

    graph = R.drawPDF()
    view = viewer.View(graph)




.. image:: /auto_reliability_sensitivity/reliability/images/sphx_glr_plot_axial_stressed_beam_quickstart_001.png
    :alt: plot axial stressed beam quickstart
    :class: sphx-glr-single-img





Our second marginal is a `Normal` univariate distribution :


.. code-block:: default

    F = sm.distribution_F








We draw the PDF of the second marginal.


.. code-block:: default

    graph = F.drawPDF()
    view = viewer.View(graph)




.. image:: /auto_reliability_sensitivity/reliability/images/sphx_glr_plot_axial_stressed_beam_quickstart_002.png
    :alt: plot axial stressed beam quickstart
    :class: sphx-glr-single-img





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 :


.. code-block:: default

    myDistribution = sm.distribution








We create a `RandomVector` from the `Distribution`, which will then be fed to the limit state function.


.. code-block:: default

    inputRandomVector = ot.RandomVector(myDistribution)








Finally we create a `CompositeRandomVector` by associating the limit state function with the input random vector.


.. code-block:: default

    outputRandomVector = ot.CompositeRandomVector(limitStateFunction, inputRandomVector)








Exact computation
-----------------

The simplest method is to perform an exact computation based on the arithmetic of distributions.


.. code-block:: default

    D = 0.02









.. code-block:: default

    G = R-F/(D**2/4 * np.pi)









.. code-block:: default

    G.computeCDF(0.)





.. rst-class:: sphx-glr-script-out

 Out:

 .. code-block:: none


    0.029198194624830504



This is the exact result from the description of this example.

Distribution of the output
--------------------------

Plot the distribution of the output.


.. code-block:: default

    sampleSize = 500
    sampleG = outputRandomVector.getSample(sampleSize)
    graph = ot.HistogramFactory().build(sampleG).drawPDF()
    view = viewer.View(graph)




.. image:: /auto_reliability_sensitivity/reliability/images/sphx_glr_plot_axial_stressed_beam_quickstart_003.png
    :alt: y0 PDF
    :class: sphx-glr-single-img





Estimate the probability with Monte-Carlo
-----------------------------------------

We first create a `ThresholdEvent` based on the output `RandomVector`, the operator and the threshold.


.. code-block:: default

    myEvent = ot.ThresholdEvent(outputRandomVector, ot.Less(), 0.0)








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`. 


.. 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)








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).


.. code-block:: default

    initialNumberOfCall = limitStateFunction.getEvaluationCallsNumber()








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.


.. code-block:: default

    algoMC.run()








We can then get the results of the algorithm. 


.. 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 = 12692
    Pf =  0.030570438071226093
    CV = 0.04998529582572751




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.


.. code-block:: default

    graph = algoMC.drawProbabilityConvergence()
    graph.setLogScale(ot.GraphImplementation.LOGX)
    view = viewer.View(graph)




.. image:: /auto_reliability_sensitivity/reliability/images/sphx_glr_plot_axial_stressed_beam_quickstart_004.png
    :alt: ProbabilitySimulationAlgorithm convergence graph at level 0.95
    :class: sphx-glr-single-img





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.


.. code-block:: default

    alpha = 0.05









.. 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.027575,0.033565]




This interval is consistent with the exact probability :math:`P_f=0.02920`.

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.455 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 <plot_axial_stressed_beam_quickstart.py>`



  .. container:: sphx-glr-download sphx-glr-download-jupyter

     :download:`Download Jupyter notebook: plot_axial_stressed_beam_quickstart.ipynb <plot_axial_stressed_beam_quickstart.ipynb>`


.. only:: html

 .. rst-class:: sphx-glr-signature

    `Gallery generated by Sphinx-Gallery <https://sphinx-gallery.github.io>`_