# Estimate a buckling probability¶

In this example, we estimate the probability that the output of a function exceeds a given threshold with the FORM method, the SORM method and an advanced sampling method.

We consider the stiffened panel model.

## Define the model¶

```from openturns.usecases import stiffened_panel
import openturns as ot
import openturns.viewer as viewer

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

We load the stiffened panel model from the usecases module :

```panel = stiffened_panel.StiffenedPanel()
distribution = panel.distribution
model = panel.model
```

See the input distribution

```distribution
```
JointDistribution
• name=JointDistribution
• dimension: 10
• description=[E (Pa),nu (-),h_c (m),ell (m),f_1 (m),f_2 (m),t (m),a (m),b_0 (m),p (m)]#10
• copula: IndependentCopula(dimension = 10)
Index Variable Distribution
0 E (Pa) TruncatedNormal(mu = 1.1e+11, sigma = 5.5e+10, a = 9.9e+10, b = 1.21e+11)
1 nu (-) Uniform(a = 0.3675, b = 0.3825)
2 h_c (m) Uniform(a = 0.0285, b = 0.0315)
3 ell (m) Uniform(a = 0.04655, b = 0.05145)
4 f_1 (m) Uniform(a = 0.0266, b = 0.0294)
5 f_2 (m) Uniform(a = 0.00627, b = 0.00693)
6 t (m) Uniform(a = 8.02e-05, b = 8.181e-05)
7 a (m) Uniform(a = 0.6039, b = 0.6161)
8 b_0 (m) Uniform(a = 0.04455, b = 0.04545)
9 p (m) Uniform(a = 0.03762, b = 0.03838)

See the model

```model.getOutputDescription()
```
[(N_{xy})_{cr}]

Draw the distribution of a sample of the output.

```sampleSize = 1000
inputSample = distribution.getSample(sampleSize)
outputSample = model(inputSample)
graph = ot.HistogramFactory().build(outputSample).drawPDF()
_ = viewer.View(graph)
```

## Define the event¶

Then we create the event whose probability we want to estimate.

```vect = ot.RandomVector(distribution)
event.setName("buckling")
```

## Estimate the probability with FORM¶

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.

```startingPoint = distribution.getMean()
algo = ot.FORM(optimAlgo, event, startingPoint)
n0 = model.getCallsNumber()
algo.run()
n1 = model.getCallsNumber()
result = algo.getResult()
standardSpaceDesignPoint = result.getStandardSpaceDesignPoint()
```

Retrieve results.

```result = algo.getResult()
probability = result.getEventProbability()
print("Pf (FORM)=%.3e" % probability, "nb evals=", n1 - n0)
```
```Pf (FORM)=7.620e-04 nb evals= 1
```

Importance factors.

```graph = result.drawImportanceFactors()
view = viewer.View(graph)
```

## Estimate the probability with SORM¶

Run SORM.

```algo = ot.SORM(optimAlgo, event, startingPoint)
n0 = model.getCallsNumber()
algo.run()
n1 = model.getCallsNumber()
```

Retrieve results.

```result = algo.getResult()
probability = result.getEventProbabilityBreitung()
print("Pf (SORM)=%.3e" % probability, "nb evals=", n1 - n0)
```
```Pf (SORM)=7.027e-05 nb evals= 1
```

We see that the FORM and SORM approximations give significantly different results. Use a simulation algorithm to get a confidence interval.

## Estimate the probability with PostAnalyticalControlledImportanceSampling¶

```algo = ot.PostAnalyticalControlledImportanceSampling(result)
algo.setBlockSize(100)
algo.setMaximumOuterSampling(100)
algo.setMaximumCoefficientOfVariation(0.1)
n0 = model.getCallsNumber()
algo.run()
n1 = model.getCallsNumber()
result = algo.getResult()
Pf = result.getProbabilityEstimate()
print("Pf (sim) = %.3e" % Pf, "nb evals=", n1 - n0)
width = result.getConfidenceLength(0.95)
print("C.I (95%)=[" + "%.3e" % (Pf - 0.5 * width), ",%.3e" % (Pf + 0.5 * width), "]")
```
```Pf (sim) = 5.914e-05 nb evals= 0
C.I (95%)=[3.060e-05 ,8.769e-05 ]
```