.. 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.py"
.. LINE NUMBERS ARE GIVEN BELOW.

.. only:: html

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

        Click :ref:`here <sphx_glr_download_auto_reliability_sensitivity_reliability_plot_axial_stressed_beam.py>`
        to download the full example code

.. rst-class:: sphx-glr-example-title

.. _sphx_glr_auto_reliability_sensitivity_reliability_plot_axial_stressed_beam.py:


Axial stressed beam : comparing different methods to estimate a probability
===========================================================================

.. GENERATED FROM PYTHON SOURCE LINES 6-13

In this example, we compare four methods to estimate the probability in the :ref:`axial stressed beam <use-case-stressed-beam>` example :

* Monte-Carlo simulation,
* FORM,
* directional sampling,
* importance sampling with FORM design point: FORM-IS.


.. GENERATED FROM PYTHON SOURCE LINES 15-17

Define the model
----------------

.. GENERATED FROM PYTHON SOURCE LINES 19-26

.. code-block:: default

    from __future__ import print_function

    import openturns as ot
    import openturns.viewer as viewer
    from matplotlib import pylab as plt
    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-31

.. code-block:: default

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








.. GENERATED FROM PYTHON SOURCE LINES 32-33

The limit state function is defined in the `model` field of the data class :

.. GENERATED FROM PYTHON SOURCE LINES 33-35

.. code-block:: default

    limitStateFunction = sm.model








.. GENERATED FROM PYTHON SOURCE LINES 36-38

The probabilistic model of the axial stressed beam is defined in the data class.
We get the first marginal and draw it :

.. GENERATED FROM PYTHON SOURCE LINES 38-42

.. code-block:: default

    R_dist = sm.distribution_R
    graph = R_dist.drawPDF()
    view = viewer.View(graph)




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





.. GENERATED FROM PYTHON SOURCE LINES 43-44

We get the second marginal and draw it :

.. GENERATED FROM PYTHON SOURCE LINES 46-50

.. code-block:: default

    F_dist = sm.distribution_F
    graph = F_dist.drawPDF()
    view = viewer.View(graph)




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





.. GENERATED FROM PYTHON SOURCE LINES 51-52

These independent marginals define the joint distribution of the input parameters :

.. GENERATED FROM PYTHON SOURCE LINES 52-55

.. code-block:: default

    myDistribution = sm.distribution









.. GENERATED FROM PYTHON SOURCE LINES 56-57

We create a `RandomVector` from the `Distribution`, then a composite random vector. Finally, we create a `ThresholdEvent` from this `RandomVector`.

.. GENERATED FROM PYTHON SOURCE LINES 59-63

.. code-block:: default

    inputRandomVector = ot.RandomVector(myDistribution)
    outputRandomVector = ot.CompositeRandomVector(limitStateFunction, inputRandomVector)
    myEvent = ot.ThresholdEvent(outputRandomVector, ot.Less(), 0.0)








.. GENERATED FROM PYTHON SOURCE LINES 64-66

Using Monte Carlo simulations
-----------------------------

.. GENERATED FROM PYTHON SOURCE LINES 68-77

.. code-block:: default

    cv = 0.05
    NbSim = 100000

    experiment = ot.MonteCarloExperiment()
    algoMC = ot.ProbabilitySimulationAlgorithm(myEvent, experiment)
    algoMC.setMaximumOuterSampling(NbSim)
    algoMC.setBlockSize(1)
    algoMC.setMaximumCoefficientOfVariation(cv)








.. GENERATED FROM PYTHON SOURCE LINES 78-79

For statistics about the algorithm

.. GENERATED FROM PYTHON SOURCE LINES 79-81

.. code-block:: default

    initialNumberOfCall = limitStateFunction.getEvaluationCallsNumber()








.. GENERATED FROM PYTHON SOURCE LINES 82-83

Perform the analysis.

.. GENERATED FROM PYTHON SOURCE LINES 85-87

.. code-block:: default

    algoMC.run()








.. GENERATED FROM PYTHON SOURCE LINES 88-95

.. code-block:: default

    result = algoMC.getResult()
    probabilityMonteCarlo = result.getProbabilityEstimate()
    numberOfFunctionEvaluationsMonteCarlo = limitStateFunction.getEvaluationCallsNumber() - initialNumberOfCall
    print('Number of calls to the limit state =', numberOfFunctionEvaluationsMonteCarlo)
    print('Pf = ', probabilityMonteCarlo)
    print('CV =', result.getCoefficientOfVariation())





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

 Out:

 .. code-block:: none

    Number of calls to the limit state = 13289
    Pf =  0.029272330498908847
    CV = 0.049954418027493265




.. GENERATED FROM PYTHON SOURCE LINES 96-100

.. 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_003.png
    :alt: ProbabilitySimulationAlgorithm convergence graph at level 0.95
    :class: sphx-glr-single-img





.. GENERATED FROM PYTHON SOURCE LINES 101-103

Using FORM analysis
-------------------

.. GENERATED FROM PYTHON SOURCE LINES 105-106

We create a NearestPoint algorithm

.. GENERATED FROM PYTHON SOURCE LINES 106-115

.. code-block:: default

    myCobyla = ot.Cobyla()
    # Resolution options:
    eps = 1e-3
    myCobyla.setMaximumEvaluationNumber(100)
    myCobyla.setMaximumAbsoluteError(eps)
    myCobyla.setMaximumRelativeError(eps)
    myCobyla.setMaximumResidualError(eps)
    myCobyla.setMaximumConstraintError(eps)








.. GENERATED FROM PYTHON SOURCE LINES 116-117

For statistics about the algorithm

.. GENERATED FROM PYTHON SOURCE LINES 117-119

.. code-block:: default

    initialNumberOfCall = limitStateFunction.getEvaluationCallsNumber()








.. GENERATED FROM PYTHON SOURCE LINES 120-121

We create a FORM algorithm. The first parameter is a NearestPointAlgorithm. The second parameter is an event. The third parameter is a starting point for the design point research.

.. GENERATED FROM PYTHON SOURCE LINES 123-125

.. code-block:: default

    algoFORM = ot.FORM(myCobyla, myEvent, myDistribution.getMean())








.. GENERATED FROM PYTHON SOURCE LINES 126-127

Perform the analysis.

.. GENERATED FROM PYTHON SOURCE LINES 129-131

.. code-block:: default

    algoFORM.run()








.. GENERATED FROM PYTHON SOURCE LINES 132-138

.. code-block:: default

    resultFORM = algoFORM.getResult()
    numberOfFunctionEvaluationsFORM = limitStateFunction.getEvaluationCallsNumber() - initialNumberOfCall
    probabilityFORM = resultFORM.getEventProbability()
    print('Number of calls to the limit state =', numberOfFunctionEvaluationsFORM)
    print('Pf =', probabilityFORM)





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

 Out:

 .. code-block:: none

    Number of calls to the limit state = 98
    Pf = 0.0299827855823147




.. GENERATED FROM PYTHON SOURCE LINES 139-142

.. code-block:: default

    graph = resultFORM.drawImportanceFactors()
    view = viewer.View(graph)




.. image:: /auto_reliability_sensitivity/reliability/images/sphx_glr_plot_axial_stressed_beam_004.png
    :alt: Importance Factors from Design Point - Unnamed
    :class: sphx-glr-single-img





.. GENERATED FROM PYTHON SOURCE LINES 143-145

Using Directional sampling
--------------------------

.. GENERATED FROM PYTHON SOURCE LINES 147-148

Resolution options:

.. GENERATED FROM PYTHON SOURCE LINES 148-156

.. code-block:: default

    cv = 0.05
    NbSim = 10000

    algoDS = ot.DirectionalSampling(myEvent)
    algoDS.setMaximumOuterSampling(NbSim)
    algoDS.setBlockSize(1)
    algoDS.setMaximumCoefficientOfVariation(cv)








.. GENERATED FROM PYTHON SOURCE LINES 157-158

For statistics about the algorithm

.. GENERATED FROM PYTHON SOURCE LINES 158-160

.. code-block:: default

    initialNumberOfCall = limitStateFunction.getEvaluationCallsNumber()








.. GENERATED FROM PYTHON SOURCE LINES 161-162

Perform the analysis.

.. GENERATED FROM PYTHON SOURCE LINES 164-166

.. code-block:: default

    algoDS.run()








.. GENERATED FROM PYTHON SOURCE LINES 167-174

.. code-block:: default

    result = algoDS.getResult()
    probabilityDirectionalSampling = result.getProbabilityEstimate()
    numberOfFunctionEvaluationsDirectionalSampling = limitStateFunction.getEvaluationCallsNumber() - initialNumberOfCall
    print('Number of calls to the limit state =', numberOfFunctionEvaluationsDirectionalSampling)
    print('Pf = ', probabilityDirectionalSampling)
    print('CV =', result.getCoefficientOfVariation())





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

 Out:

 .. code-block:: none

    Number of calls to the limit state = 10139
    Pf =  0.0276771208465632
    CV = 0.04988350600378961




.. GENERATED FROM PYTHON SOURCE LINES 175-179

.. code-block:: default

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




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





.. GENERATED FROM PYTHON SOURCE LINES 180-182

Using importance sampling with FORM design point: FORM-IS
---------------------------------------------------------

.. GENERATED FROM PYTHON SOURCE LINES 184-185

The `getStandardSpaceDesignPoint` method returns the design point in the U-space.

.. GENERATED FROM PYTHON SOURCE LINES 187-190

.. code-block:: default

    standardSpaceDesignPoint = resultFORM.getStandardSpaceDesignPoint()
    standardSpaceDesignPoint






.. raw:: html

    <div class="output_subarea output_html rendered_html output_result">
    <p>[-1.59355,0.999463]</p>
    </div>
    <br />
    <br />

.. GENERATED FROM PYTHON SOURCE LINES 191-192

The key point is to define the importance distribution in the U-space. To define it, we use a multivariate standard Gaussian and configure it so that the center is equal to the design point in the U-space.

.. GENERATED FROM PYTHON SOURCE LINES 194-197

.. code-block:: default

    dimension = myDistribution.getDimension()
    dimension





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

 Out:

 .. code-block:: none


    2



.. GENERATED FROM PYTHON SOURCE LINES 198-202

.. code-block:: default

    myImportance = ot.Normal(dimension)
    myImportance.setMean(standardSpaceDesignPoint)
    myImportance






.. raw:: html

    <div class="output_subarea output_html rendered_html output_result">
    <p>Normal(mu = [-1.59355,0.999463], sigma = [1,1], R = [[ 1 0 ]<br>
     [ 0 1 ]])</p>
    </div>
    <br />
    <br />

.. GENERATED FROM PYTHON SOURCE LINES 203-204

Create the design of experiment corresponding to importance sampling. This generates a `WeightedExperiment` with weights corresponding to the importance distribution.

.. GENERATED FROM PYTHON SOURCE LINES 206-208

.. code-block:: default

    experiment = ot.ImportanceSamplingExperiment(myImportance)








.. GENERATED FROM PYTHON SOURCE LINES 209-210

Create the standard event corresponding to the event. This transforms the original problem into the U-space, with Gaussian independent marginals.

.. GENERATED FROM PYTHON SOURCE LINES 212-214

.. code-block:: default

    standardEvent = ot.StandardEvent(myEvent)








.. GENERATED FROM PYTHON SOURCE LINES 215-216

We then create the simulation algorithm. 

.. GENERATED FROM PYTHON SOURCE LINES 218-222

.. code-block:: default

    algo = ot.ProbabilitySimulationAlgorithm(standardEvent, experiment)
    algo.setMaximumCoefficientOfVariation(cv)
    algo.setMaximumOuterSampling(40000)








.. GENERATED FROM PYTHON SOURCE LINES 223-224

For statistics about the algorithm

.. GENERATED FROM PYTHON SOURCE LINES 224-226

.. code-block:: default

    initialNumberOfCall = limitStateFunction.getEvaluationCallsNumber()








.. GENERATED FROM PYTHON SOURCE LINES 227-229

.. code-block:: default

    algo.run()








.. GENERATED FROM PYTHON SOURCE LINES 230-231

retrieve results

.. GENERATED FROM PYTHON SOURCE LINES 231-238

.. code-block:: default

    result = algo.getResult()
    probabilityFORMIS = result.getProbabilityEstimate()
    numberOfFunctionEvaluationsFORMIS = limitStateFunction.getEvaluationCallsNumber() - initialNumberOfCall
    print('Number of calls to the limit state =', numberOfFunctionEvaluationsFORMIS)
    print('Pf = ', probabilityFORMIS)
    print('CV =', result.getCoefficientOfVariation())





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

 Out:

 .. code-block:: none

    Number of calls to the limit state = 883
    Pf =  0.02985232966804074
    CV = 0.04991286447957941




.. GENERATED FROM PYTHON SOURCE LINES 239-241

Conclusion
----------

.. GENERATED FROM PYTHON SOURCE LINES 243-244

We now compare the different methods in terms of accuracy and speed.

.. GENERATED FROM PYTHON SOURCE LINES 246-249

.. code-block:: default

    import numpy as np









.. GENERATED FROM PYTHON SOURCE LINES 250-251

The following function computes the number of correct base-10 digits in the computed result compared to the exact result.

.. GENERATED FROM PYTHON SOURCE LINES 253-258

.. code-block:: default

    def computeLogRelativeError(exact, computed):
        logRelativeError = -np.log10(abs(exact - computed) / abs(exact))
        return logRelativeError









.. GENERATED FROM PYTHON SOURCE LINES 259-260

The following function prints the results.

.. GENERATED FROM PYTHON SOURCE LINES 262-275

.. code-block:: default

    def printMethodSummary(name, computedProbability, numberOfFunctionEvaluations):
        print("---")
        print(name,":")
        print('Number of calls to the limit state =', numberOfFunctionEvaluations)
        print('Pf = ', computedProbability)
        exactProbability = 0.02919819462483051
        logRelativeError = computeLogRelativeError(exactProbability, computedProbability)
        print("Number of correct digits=%.3f" % (logRelativeError))
        performance = logRelativeError/numberOfFunctionEvaluations
        print("Performance=%.2e (correct digits/evaluation)" % (performance))
        return









.. GENERATED FROM PYTHON SOURCE LINES 276-281

.. code-block:: default

    printMethodSummary("Monte-Carlo", probabilityMonteCarlo, numberOfFunctionEvaluationsMonteCarlo)
    printMethodSummary("FORM", probabilityFORM, numberOfFunctionEvaluationsFORM)
    printMethodSummary("DirectionalSampling", probabilityDirectionalSampling, numberOfFunctionEvaluationsDirectionalSampling)
    printMethodSummary("FORM-IS", probabilityFORMIS, numberOfFunctionEvaluationsFORMIS)





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

 Out:

 .. code-block:: none

    ---
    Monte-Carlo :
    Number of calls to the limit state = 13289
    Pf =  0.029272330498908847
    Number of correct digits=2.595
    Performance=1.95e-04 (correct digits/evaluation)
    ---
    FORM :
    Number of calls to the limit state = 98
    Pf =  0.0299827855823147
    Number of correct digits=1.571
    Performance=1.60e-02 (correct digits/evaluation)
    ---
    DirectionalSampling :
    Number of calls to the limit state = 10139
    Pf =  0.0276771208465632
    Number of correct digits=1.283
    Performance=1.27e-04 (correct digits/evaluation)
    ---
    FORM-IS :
    Number of calls to the limit state = 883
    Pf =  0.02985232966804074
    Number of correct digits=1.650
    Performance=1.87e-03 (correct digits/evaluation)




.. GENERATED FROM PYTHON SOURCE LINES 282-288

We see that all three methods produce the correct probability, but not with the same accuracy. In this case, we have found the correct order of magnitude of the probability, i.e. between one and two correct digits. There is, however, a significant difference in computational performance (measured here by the number of function evaluations).

* The fastest method is the FORM method, which produces more than 1 correct digit with less than 98 function evaluations with a performance equal to :math:`1.60 \times 10^{-2}` (correct digits/evaluation). A practical limitation is that the FORM method does not produce a confidence interval: there is no guarantee that the computed probability is correct.
* The slowest method is Monte-Carlo simulation, which produces more than 1 correct digit with 12806 function evaluations. This is associated with a very slow performance equal to :math:`1.11 \times 10^{-4}` (correct digits/evaluation). The interesting point with the Monte-Carlo simulation is that the method produces a confidence interval.
* The DirectionalSampling method is somewhat in-between the two previous methods.
* The FORM-IS method produces 2 correct digits and has a small number of function evaluations. It has an intermediate performance equal to :math:`2.37\times 10^{-3}` (correct digits/evaluation). It combines the best of the both worlds: it has the small number of function evaluation of FORM computation and the confidence interval of Monte-Carlo simulation.


.. rst-class:: sphx-glr-timing

   **Total running time of the script:** ( 0 minutes  0.852 seconds)


.. _sphx_glr_download_auto_reliability_sensitivity_reliability_plot_axial_stressed_beam.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.py <plot_axial_stressed_beam.py>`



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

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


.. only:: html

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

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