# Estimate a probability with Monte Carlo¶

In this example we estimate a probability by means of a simulation algorithm, the Monte-Carlo algorithm. To do this, we need the classes MonteCarloExperiment and ProbabilitySimulationAlgorithm.

## Introduction¶

We consider a simple beam stressed by a traction load F at both sides.

The geometry is supposed to be deterministic; the diameter D is equal to:

By definition, the yield stress is the load divided by the surface. Since the surface is , the stress is:

Failure occurs when the beam plastifies, i.e. when the axial stress gets larger than the yield stress:

where is the strength.

Therefore, the limit state function is:

for any .

The value of the parameter is such that:

We consider the following distribution functions.

Variable

Distribution

R

LogNormal(, ) [Pa]

F

Normal(, ) [N]

where and are the mean and the variance of .

The failure probability is:

The exact is

[1]:

from __future__ import print_function
import openturns as ot


Create the joint distribution of the parameters.

[2]:

distribution_R = ot.LogNormalMuSigma(300.0, 30.0, 0.0).getDistribution()
distribution_F = ot.Normal(75e3, 5e3)
marginals = [distribution_R, distribution_F]
distribution = ot.ComposedDistribution(marginals)


Create the model.

[3]:

model = ot.SymbolicFunction(['R', 'F'], ['R-F/(pi_*100.0)'])


Create the event whose probability we want to estimate.

[4]:

vect = ot.RandomVector(distribution)
G = ot.CompositeRandomVector(model, vect)
event = ot.ThresholdEvent(G, ot.Less(), 0.0)


Create a Monte Carlo algorithm.

[5]:

experiment = ot.MonteCarloExperiment()
algo = ot.ProbabilitySimulationAlgorithm(event, experiment)
algo.setMaximumCoefficientOfVariation(0.05)
algo.setMaximumOuterSampling(int(1e5))
algo.run()


Retrieve results.

[6]:

result = algo.getResult()
probability = result.getProbabilityEstimate()
print('Pf=', probability)

Pf= 0.03029829767296579