Exploitation of simulation algorithm resultsΒΆ

In this example we are going to retrieve all the results proposed by probability simulation algorithms:

  • the probability estimate

  • the estimator variance

  • the confidence interval

  • the convergence graph of the estimator

  • the stored input and output numerical samples

  • importance factors

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

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

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)"])
modelCallNumberBefore = model.getEvaluationCallsNumber()
modelGradientCallNumberBefore = model.getGradientCallsNumber()
modelHessianCallNumberBefore = model.getHessianCallsNumber()

To have access to the input and output samples after the simulation, activate the History mechanism.

model = ot.MemoizeFunction(model)

Remove all the values stored in the history mechanism. Care : it is done regardless the status of the History mechanism.

model.clearHistory()

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)
algo.setMaximumCoefficientOfVariation(0.1)
algo.setMaximumStandardDeviation(0.001)
algo.setMaximumOuterSampling(int(1e4))

Define the HistoryStrategy to store the values of P_n and \sigma_n used ot draw the convergence graph. Compact strategy : N points

N = 1000
algo.setConvergenceStrategy(ot.Compact(N))
algo.run()

Retrieve result structure.

result = algo.getResult()

Display the simulation event probability.

result.getProbabilityEstimate()
0.03184504267892307

Criteria 3 : Display the Standard Deviation of the estimator

result.getStandardDeviation()
0.0031814741942203543

Display the variance of the simulation probability estimator.

result.getVarianceEstimate()
1.0121778048490053e-05

Criteria 2 : Display the number of iterations of the simulation

result.getOuterSampling()
3046

Display the total number of evaluations of the model

result.getOuterSampling() * result.getBlockSize()
3046

Save the number of calls to the model, its gradient and hessian done so far.

modelCallNumberAfter = model.getEvaluationCallsNumber()
modelGradientCallNumberAfter = model.getGradientCallsNumber()
modelHessianCallNumberAfter = model.getHessianCallsNumber()

Display the number of iterations executed and the number of evaluations of the model.

modelCallNumberAfter - modelCallNumberBefore
3046

Get the mean point in event domain care : only for Monte Carlo and LHS sampling methods.

result.getMeanPointInEventDomain()

[247.822,80985.1]



Get the associated importance factors care : only for Monte Carlo and LHS sampling methods.

result.getImportanceFactors()

[X0 : 0.708375, X1 : 0.291625]



graph = result.drawImportanceFactors()
view = viewer.View(graph)
Importance Factors from Simulation - v0

Display the confidence interval length centered around the MonteCarlo probability. The confidence interval is

IC = [\tilde{p} - 0.5 \ell, \tilde{p} + 0.5 \ell]

with level 0.95, where \tilde{p} is the estimated probability and \ell is the confidence interval length.

probability = result.getProbabilityEstimate()
length95 = result.getConfidenceLength(0.95)
print("0.95 Confidence Interval length = ", length95)
print(
    "IC at 0.95 = [",
    probability - 0.5 * length95,
    "; ",
    probability + 0.5 * length95,
    "]",
)
0.95 Confidence Interval length =  0.012471149676830965
IC at 0.95 = [ 0.025609467840507584 ;  0.03808061751733855 ]

Draw the convergence graph and the confidence interval of level alpha. By default, alpha = 0.95.

alpha = 0.90
graph = algo.drawProbabilityConvergence(alpha)
view = viewer.View(graph)
ProbabilitySimulationAlgorithm convergence graph at level 0.9

Get the numerical samples of the input and output random vectors stored according to the History Strategy specified and used to evaluate the probability estimator and its variance.

inputSampleStored = model.getInputHistory()
outputSampleStored = model.getOutputHistory()
inputSampleStored
v0v1
0290.454686759.83
1315.177975776.36
2348.053777994.47
...
3043322.057271236.52
3044294.31877531.99
3045262.282385847.2


Get the values of the estimator and its variance stored according to the History Strategy specified and used to draw the convergence graph.

estimator_probability_sample = algo.getConvergenceStrategy().getSample()[0]
estimator_variance_sample = algo.getConvergenceStrategy().getSample()[1]
print(estimator_probability_sample, estimator_variance_sample)
plt.show()
[0,-1] [0,-1]