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

        :ref:`Go to the end <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 7-14

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 16-18

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

.. GENERATED FROM PYTHON SOURCE LINES 20-28

.. code-block:: Python

    import numpy as np
    from openturns.usecases import stressed_beam

    import openturns as ot
    import openturns.viewer as viewer

    ot.Log.Show(ot.Log.NONE)








.. GENERATED FROM PYTHON SOURCE LINES 29-30

We load the model from the usecases module :

.. GENERATED FROM PYTHON SOURCE LINES 30-32

.. code-block:: Python

    sm = stressed_beam.AxialStressedBeam()








.. GENERATED FROM PYTHON SOURCE LINES 33-34

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

.. GENERATED FROM PYTHON SOURCE LINES 34-36

.. code-block:: Python

    limitStateFunction = sm.model








.. GENERATED FROM PYTHON SOURCE LINES 37-39

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 39-43

.. code-block:: Python

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




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





.. GENERATED FROM PYTHON SOURCE LINES 44-45

We get the second marginal and draw it :

.. GENERATED FROM PYTHON SOURCE LINES 47-51

.. code-block:: Python

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




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





.. GENERATED FROM PYTHON SOURCE LINES 52-53

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

.. GENERATED FROM PYTHON SOURCE LINES 53-56

.. code-block:: Python

    myDistribution = sm.distribution









.. GENERATED FROM PYTHON SOURCE LINES 57-58

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 60-64

.. code-block:: Python

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








.. GENERATED FROM PYTHON SOURCE LINES 65-67

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

.. GENERATED FROM PYTHON SOURCE LINES 69-78

.. code-block:: Python

    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 79-80

For statistics about the algorithm

.. GENERATED FROM PYTHON SOURCE LINES 80-82

.. code-block:: Python

    initialNumberOfCall = limitStateFunction.getEvaluationCallsNumber()








.. GENERATED FROM PYTHON SOURCE LINES 83-84

Perform the analysis.

.. GENERATED FROM PYTHON SOURCE LINES 86-88

.. code-block:: Python

    algoMC.run()








.. GENERATED FROM PYTHON SOURCE LINES 89-98

.. code-block:: Python

    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

 .. code-block:: none

    Number of calls to the limit state = 12927
    Pf =  0.03001469791908405
    CV = 0.049999621186217445




.. GENERATED FROM PYTHON SOURCE LINES 99-103

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





.. GENERATED FROM PYTHON SOURCE LINES 104-106

Using LHS simulation
--------------------

.. GENERATED FROM PYTHON SOURCE LINES 106-113

.. code-block:: Python

    experiment = ot.LHSExperiment()
    experiment.setAlwaysShuffle(True)
    algo = ot.ProbabilitySimulationAlgorithm(myEvent, experiment)
    algo.setMaximumOuterSampling(NbSim)
    algo.setBlockSize(1)
    algo.setMaximumCoefficientOfVariation(cv)








.. GENERATED FROM PYTHON SOURCE LINES 114-115

For statistics about the algorithm

.. GENERATED FROM PYTHON SOURCE LINES 115-117

.. code-block:: Python

    initialNumberOfCall = limitStateFunction.getEvaluationCallsNumber()








.. GENERATED FROM PYTHON SOURCE LINES 118-119

Perform the analysis.

.. GENERATED FROM PYTHON SOURCE LINES 121-123

.. code-block:: Python

    algo.run()








.. GENERATED FROM PYTHON SOURCE LINES 124-134

.. code-block:: Python

    resultLHS = algo.getResult()

    numberOfFunctionEvaluationsLHS = (
        limitStateFunction.getEvaluationCallsNumber() - initialNumberOfCall
    )
    probabilityLHS = result.getProbabilityEstimate()
    print("Number of calls to the limit state =", numberOfFunctionEvaluationsLHS)
    print("Pf = ", probabilityLHS)
    print("CV =", result.getCoefficientOfVariation())





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

 .. code-block:: none

    Number of calls to the limit state = 13532
    Pf =  0.03001469791908405
    CV = 0.049999621186217445




.. GENERATED FROM PYTHON SOURCE LINES 135-140

.. code-block:: Python

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





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





.. GENERATED FROM PYTHON SOURCE LINES 141-143

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

.. GENERATED FROM PYTHON SOURCE LINES 145-146

We create a NearestPoint algorithm

.. GENERATED FROM PYTHON SOURCE LINES 146-155

.. code-block:: Python

    algoOptim = ot.AbdoRackwitz()
    # Resolution options:
    eps = 1e-3
    algoOptim.setMaximumCallsNumber(1000)
    algoOptim.setMaximumAbsoluteError(eps)
    algoOptim.setMaximumRelativeError(eps)
    algoOptim.setMaximumResidualError(eps)
    algoOptim.setMaximumConstraintError(eps)








.. GENERATED FROM PYTHON SOURCE LINES 156-157

For statistics about the algorithm

.. GENERATED FROM PYTHON SOURCE LINES 157-159

.. code-block:: Python

    initialNumberOfCall = limitStateFunction.getEvaluationCallsNumber()








.. GENERATED FROM PYTHON SOURCE LINES 160-161

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 163-165

.. code-block:: Python

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








.. GENERATED FROM PYTHON SOURCE LINES 166-167

Perform the analysis.

.. GENERATED FROM PYTHON SOURCE LINES 169-171

.. code-block:: Python

    algoFORM.run()








.. GENERATED FROM PYTHON SOURCE LINES 172-180

.. code-block:: Python

    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

 .. code-block:: none

    Number of calls to the limit state = 7
    Pf = 0.029982795362833412




.. GENERATED FROM PYTHON SOURCE LINES 181-184

.. code-block:: Python

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




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





.. GENERATED FROM PYTHON SOURCE LINES 185-187

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

.. GENERATED FROM PYTHON SOURCE LINES 189-190

Resolution options:

.. GENERATED FROM PYTHON SOURCE LINES 190-198

.. code-block:: Python

    cv = 0.05
    NbSim = 10000

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








.. GENERATED FROM PYTHON SOURCE LINES 199-200

For statistics about the algorithm

.. GENERATED FROM PYTHON SOURCE LINES 200-202

.. code-block:: Python

    initialNumberOfCall = limitStateFunction.getEvaluationCallsNumber()








.. GENERATED FROM PYTHON SOURCE LINES 203-204

Perform the analysis.

.. GENERATED FROM PYTHON SOURCE LINES 206-208

.. code-block:: Python

    algoDS.run()








.. GENERATED FROM PYTHON SOURCE LINES 209-221

.. code-block:: Python

    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

 .. code-block:: none

    Number of calls to the limit state = 9479
    Pf =  0.02874805526680341
    CV = 0.04988071645434962




.. GENERATED FROM PYTHON SOURCE LINES 222-226

.. code-block:: Python

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




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





.. GENERATED FROM PYTHON SOURCE LINES 227-229

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

.. GENERATED FROM PYTHON SOURCE LINES 231-232

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

.. GENERATED FROM PYTHON SOURCE LINES 234-237

.. code-block:: Python

    standardSpaceDesignPoint = resultFORM.getStandardSpaceDesignPoint()
    standardSpaceDesignPoint






.. raw:: html

    <div class="output_subarea output_html rendered_html output_result">
    class=Point name=Standard Space Design Point dimension=2 values=[-1.59396,0.99882]
    </div>
    <br />
    <br />

.. GENERATED FROM PYTHON SOURCE LINES 238-240

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 242-245

.. code-block:: Python

    dimension = myDistribution.getDimension()
    dimension





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

 .. code-block:: none


    2



.. GENERATED FROM PYTHON SOURCE LINES 246-250

.. code-block:: Python

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






.. raw:: html

    <div class="output_subarea output_html rendered_html output_result">
    Normal
    <ul>
      <li>name=Normal</li>
      <li>dimension=2</li>
      <li>weight=1</li>
      <li>range=]-inf (-9.24458), (6.05667) +inf[
    ]-inf (-6.65181), (8.64945) +inf[</li>
      <li>description=[X0,X1]</li>
      <li>isParallel=true</li>
      <li>isCopula=false</li>
    </ul>

    </div>
    <br />
    <br />

.. GENERATED FROM PYTHON SOURCE LINES 251-252

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 254-256

.. code-block:: Python

    experiment = ot.ImportanceSamplingExperiment(myImportance)








.. GENERATED FROM PYTHON SOURCE LINES 257-258

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 260-262

.. code-block:: Python

    standardEvent = ot.StandardEvent(myEvent)








.. GENERATED FROM PYTHON SOURCE LINES 263-264

We then create the simulation algorithm.

.. GENERATED FROM PYTHON SOURCE LINES 266-270

.. code-block:: Python

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








.. GENERATED FROM PYTHON SOURCE LINES 271-272

For statistics about the algorithm

.. GENERATED FROM PYTHON SOURCE LINES 272-274

.. code-block:: Python

    initialNumberOfCall = limitStateFunction.getEvaluationCallsNumber()








.. GENERATED FROM PYTHON SOURCE LINES 275-277

.. code-block:: Python

    algo.run()








.. GENERATED FROM PYTHON SOURCE LINES 278-279

retrieve results

.. GENERATED FROM PYTHON SOURCE LINES 279-288

.. code-block:: Python

    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

 .. code-block:: none

    Number of calls to the limit state = 891
    Pf =  0.030668150809484762
    CV = 0.04999634278115516




.. GENERATED FROM PYTHON SOURCE LINES 289-291

Conclusion
----------

.. GENERATED FROM PYTHON SOURCE LINES 293-294

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

.. GENERATED FROM PYTHON SOURCE LINES 299-300

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 303-308

.. code-block:: Python

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









.. GENERATED FROM PYTHON SOURCE LINES 309-310

The following function prints the results.

.. GENERATED FROM PYTHON SOURCE LINES 313-326

.. code-block:: Python

    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 327-339

.. code-block:: Python

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





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

 .. code-block:: none

    ---
    Monte-Carlo :
    Number of calls to the limit state = 12927
    Pf =  0.03001469791908405
    Number of correct digits=1.553
    Performance=1.20e-04 (correct digits/evaluation)
    ---
    LHS :
    Number of calls to the limit state = 13532
    Pf =  0.03001469791908405
    Number of correct digits=1.553
    Performance=1.15e-04 (correct digits/evaluation)
    ---
    FORM :
    Number of calls to the limit state = 7
    Pf =  0.029982795362833412
    Number of correct digits=1.571
    Performance=2.24e-01 (correct digits/evaluation)
    ---
    DirectionalSampling :
    Number of calls to the limit state = 9479
    Pf =  0.02874805526680341
    Number of correct digits=1.812
    Performance=1.91e-04 (correct digits/evaluation)
    ---
    FORM-IS :
    Number of calls to the limit state = 891
    Pf =  0.030668150809484762
    Number of correct digits=1.298
    Performance=1.46e-03 (correct digits/evaluation)




.. GENERATED FROM PYTHON SOURCE LINES 340-354

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


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

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

      :download:`Download zipped: plot_axial_stressed_beam.zip <plot_axial_stressed_beam.zip>`