Parametrization of a simulation algorithmΒΆ

In this example we are going to parameterize a simulation algorithm:

  • parameters linked to the number of points generated

  • the precision of the probability estimator

  • the sample storage strategy

  • using callbacks to monitor progress and stopping criteria.

from __future__ import print_function
import openturns as ot
import openturns.viewer as viewer
from matplotlib import pylab as plt

Create the joint distribution of the parameters.

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.

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

Create the event whose probability we want to estimate.

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

Create a Monte Carlo algorithm.

experiment = ot.MonteCarloExperiment()
algo = ot.ProbabilitySimulationAlgorithm(event, experiment)

Criteria 1: Define the Maximum Coefficient of variation of the probability estimator.


Criteria 2: Define the number of iterations of the simulation.


The block size parameter represents the number of samples evaluated per iteration, useful for parallelization.


HistoryStrategy to store the values of the probability used to draw the convergence graph.

Null strategy


# Full strategy

# Compact strategy: N points
N = 1000

Use a callback to display the progress every 10%.

def progress(p):
    if p >= progress.t:
        progress.t += 10.0
        print('progress=', p, '%')
    return False
progress.t = 10.0

Use a callback to stop the simulation.

def stop():
    # here we never stop, but we could
    return False


progress= 10.0 %
progress= 20.0 %
progress= 30.0 %
progress= 40.0 %
progress= 50.0 %
progress= 60.0 %
progress= 70.0 %
progress= 80.0 %

Retrieve results.

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


Pf= 0.03252214870472134

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

Gallery generated by Sphinx-Gallery