Use the Importance Sampling algorithm

In this example we estimate a failure probability with the importance sampling simulation algorithm provided by the ImportanceSamplingExperiment class.

The main steps of this method are:

• run a FORM analysis,

• create an importance distribution based on the results of the FORM results,

• run a sampling-based probability estimate algorithm.

We shall consider the analytical example of a cantilever beam.

Define the cantilever beam model

from openturns.usecases import cantilever_beam
import openturns as ot
import openturns.viewer as otv

The cantilever beam example can be accessed in the usecases module :

cb = cantilever_beam.CantileverBeam()

The joint probability distribution of the input parameters is stored in the object cb :

distribution = cb.distribution

We create the model.

model = cb.model

Define the event

We create the event whose probability we want to estimate.

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

Run a FORM analysis

Define a solver

optimAlgo = ot.Cobyla()
optimAlgo.setMaximumCallsNumber(1000)
optimAlgo.setMaximumAbsoluteError(1.0e-10)
optimAlgo.setMaximumRelativeError(1.0e-10)
optimAlgo.setMaximumResidualError(1.0e-10)
optimAlgo.setMaximumConstraintError(1.0e-10)

Run FORM

optimAlgo.setStartingPoint(distribution.getMean())
algo = ot.FORM(optimAlgo, event)
algo.run()
result = algo.getResult()

Configure an importance sampling algorithm

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

standardSpaceDesignPoint = result.getStandardSpaceDesignPoint()
standardSpaceDesignPoint

class=Point name=Standard Space Design Point dimension=4 values=[-0.665643,4.31264,1.23029,-1.3689]

The key point is to define the importance distribution in the U-space. To define it, we use a multivariate standard Gaussian centered around the design point in the U-space.

dimension = distribution.getDimension()
dimension

4

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

Normal

• name=Normal
• dimension=4
• weight=1
• range=]-inf (-8.31627), (6.98498) +inf[ ]-inf (-3.33799), (11.9633) +inf[ ]-inf (-6.42034), (8.88091) +inf[ ]-inf (-9.01953), (6.28172) +inf[
• description=[X0,X1,X2,X3]
• isParallel=true
• isCopula=false

Create the design of experiment corresponding to importance sampling. This generates a WeightedExperiment with weights fitting to the importance distribution.

experiment = ot.ImportanceSamplingExperiment(myImportance)
type(experiment)

Create the standard event corresponding to the event. This pushes the original problem to the U-space, with Gaussian independent marginals.

standardEvent = ot.StandardEvent(event)

Run the importance sampling simulation

We then create the simulation algorithm.

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

We can retrieve results of this estimate :

result = algo.getResult()
probability = result.getProbabilityEstimate()
print("Probability = ", probability)

Probability =  4.554227632606641e-07

The confidence interval of level is:

ci95 = result.getProbabilityDistribution().computeBilateralConfidenceInterval(0.95)
print("Confidence interval (0.95) = ", ci95)

Confidence interval (0.95) =  [3.66496e-07, 5.4435e-07]

We can observe the convergence history of the estimate with the drawProbabilityConvergence method which displays the estimate and confidence interval evolution.

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

otv.View.ShowAll()