# Estimate a probability using randomized QMC¶

In this example we are going to estimate a failure probability.

## Introduction¶

Let us consider the following analytical example of a cantilever beam with Young’s modulus E, length L and section modulus I.

One end of the cantilever beam is built in a wall and we apply a concentrated bending load F at the other end of the beam, resulting in a deviation: Failure occurs when the beam deviation is too large: Four independent random variables are considered:

• E: Young’s modulus [Pa]

• L: length [m]

• I: section [m^4]

Stochastic model (simplified model, no units):

• E ~ Beta(0.93, 2.27, 2.8e7, 4.8e7)

• F ~ LogNormal(30000, 9000, 15000)

• L ~ Uniform(250, 260)

• I ~ Beta(2.5, 1.5, 3.1e2, 4.5e2)

:

from __future__ import print_function
import openturns as ot


Create the marginal distributions of the parameters.

:

dist_E = ot.Beta(0.93, 2.27, 2.8e7, 4.8e7)
dist_F = ot.LogNormalMuSigma(30000, 9000, 15000).getDistribution()
dist_L = ot.Uniform(250, 260)
dist_I = ot.Beta(2.5, 1.5, 3.1e2, 4.5e2)
marginals = [dist_E, dist_F, dist_L, dist_I]


Create the Copula.

:

RS = ot.CorrelationMatrix(4)
RS[2, 3] = -0.2
# Evaluate the correlation matrix of the Normal copula from RS
R = ot.NormalCopula.GetCorrelationFromSpearmanCorrelation(RS)
# Create the Normal copula parametrized by R
copula = ot.NormalCopula(R)


Create the joint probability distribution.

:

distribution = ot.ComposedDistribution(marginals, copula)


Create the model.

:

model = ot.SymbolicFunction(['E', 'F', 'L', 'I'], ['F*L^3/(3*E*I)'])


Create the event whose probability we want to estimate.

:

vect = ot.RandomVector(distribution)
G = ot.CompositeRandomVector(model, vect)
event = ot.ThresholdEvent(G, ot.Greater(), 30.0)


Define the low discrepancy sequence.

:

sequence = ot.SobolSequence()


Create a simulation algorithm.

:

experiment = ot.LowDiscrepancyExperiment(sequence, 1)
experiment.setRandomize(True)
algo = ot.ProbabilitySimulationAlgorithm(event, experiment)
algo.setMaximumCoefficientOfVariation(0.05)
algo.setMaximumOuterSampling(int(1e5))
algo.run()


Retrieve results.

:

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

Pf= 0.005667739454871699