Test the design point with the Strong Maximum Test

In this example we are going to validate a FORM estimation using the Strong Maximum Test.

The Strong Maximum Test helps to evaluate the quality of the design point resulting from the optimization algorithm. It checks whether the design point computed is:

• the true design point, which means a global maximum point,

• a strong design point, which means that there is no other local maximum located on the event boundary and which likelihood is slightly inferior to the design point one.

This verification is very important in order to give sense to the FORM and SORM approximations.

We briefly recall here the main principles of the test. The objective is to detect all the points in the ball of radius which are potentially the real design point (case of ) or which contribution to is not negligeable as regards the approximations Form and Sorm (case of ). The contribution of a point is considered as negligeable when its likelihood in the -space is more than -times lesser than the design point one. The radius is the distance to the -space center upon which points are considered as negligeable in the evaluation of .

from __future__ import print_function
import openturns as ot
import openturns.viewer as viewer
from matplotlib import pylab as plt
ot.Log.Show(ot.Log.NONE)

Create the model Y = x1 + 2*x2

model = ot.SymbolicFunction(["x1", "x2"], ["x1+2*x2"])

# Create the input distribution and random vector X
inputDist = ot.Normal(2)
inputDist.setDescription(['X1', 'X2'])

inputVector = ot.RandomVector(inputDist)

# Create the output random vector Y=model(X)
output = ot.CompositeRandomVector(model, inputVector)
output.setName("MyOutputY")

Create the physical event Y > 4

threshold = 4
myEvent = ot.ThresholdEvent(output, ot.Greater(), threshold)

# Create the associated standard event in the standard space
myStandardEvent = ot.StandardEvent(myEvent)

First : FORM analyses to get the design point

myCobyla = ot.Cobyla()
myStartingPoint = inputDist.getMean()
myAlgoFORM = ot.FORM(myCobyla, myEvent, myStartingPoint)
myAlgoFORM.run()
FORMResult = myAlgoFORM.getResult()
standardSpaceDesignPoint = FORMResult.getStandardSpaceDesignPoint()

Fix the importance level epsilon of the test Care : 0<epsilon<1

importanceLevel = 0.15

# Fix the accuracy level tau of the test
# Care : tau >0
# It is recommended that tau <4
accuracyLevel = 3.0

# Fix the confidence level (1-q) of the test
confidenceLevel = 0.99

# Create the Strong Maximum Test
# CARE : the event must be declared in the standard space
# 1. From the confidenceLevel parameter
mySMT_CL = ot.StrongMaximumTest(myStandardEvent, standardSpaceDesignPoint,  importanceLevel, accuracyLevel,  confidenceLevel)

# 2. Or from the  maximum number of points sampling the sphere
pointsNumber = 1000
mySMT_PN = ot.StrongMaximumTest(myStandardEvent, standardSpaceDesignPoint,  importanceLevel, accuracyLevel,  pointsNumber)

# Perform the test
mySMT_CL.run()
mySMT_PN.run()

# Get (or evaluate) the confidence level
# associated to the number of points used to sample the sphere
print('Confidence level = ', mySMT_CL.getConfidenceLevel())

# Get (or evaluate) the number of points used to sample the sphere
# associated the confidence level used
print('Points Number = ', mySMT_CL.getPointNumber())

Out:

Confidence level =  0.99
Points Number =  13

# Get all the points verifying the event  and outside the design point vicinity
# Get also the values of limit state function at these points
potentialDesignPoints = mySMT_CL.getFarDesignPointVerifyingEventPoints()
values = mySMT_CL.getFarDesignPointVerifyingEventValues()
print('Potential design points = ', potentialDesignPoints)
print('Model values = ', values)

# Get all the points verifying the event and inside the design point vicinity
# Get also the values of limit state function at these points
vicinityDesignPoint = mySMT_CL.getNearDesignPointVerifyingEventPoints()
values = mySMT_CL.getNearDesignPointVerifyingEventValues()
print('Points verifying the Event in the vicinity of the design points = ', vicinityDesignPoint)
print('Model values = ', values)

# Get all the points not verifying the event and outside the design point vicinity
# Get also the values of limit state function at these points
farSecurityPoints = mySMT_CL.getFarDesignPointViolatingEventPoints()
values = mySMT_CL.getFarDesignPointViolatingEventValues()
print('Points NOT verifying the Event outside the vicinity of the design points = ', farSecurityPoints)
print('Model values = ', values)

# Get  all the points not verifying the event and inside the design point vicinity
# Get also the values of limit state function at these points
vicinitySecurityPoints = mySMT_CL.getNearDesignPointViolatingEventPoints()
values = mySMT_CL.getNearDesignPointViolatingEventValues()
print('Points NOT verifying the Event outside the vicinity of the design points = ', vicinitySecurityPoints)
print('Model values = ', values)

Out:

Potential design points =
Model values =
Points verifying the Event in the vicinity of the design points =  0 : [ 0.790233 4.284    ]
Model values =      [ y0      ]
0 : [ 9.35824 ]
Points NOT verifying the Event outside the vicinity of the design points =   0 : [ -4.27743   0.825073 ]
 1 : [ -2.99162  -3.1666   ]
 2 : [  4.13624  -1.36699  ]
 3 : [ -4.03395  -1.64452  ]
 4 : [ -3.93906   1.86036  ]
 5 : [ -3.30804   2.83443  ]
 6 : [  3.67676  -2.33636  ]
 7 : [  2.7208   -3.40212  ]
 8 : [  1.38001  -4.13192  ]
 9 : [ -2.20529  -3.75684  ]
10 : [ -4.15565  -1.3068   ]
11 : [ -4.26232   0.899874 ]
Model values =       [ y0        ]
 0 : [ -2.62729  ]
 1 : [ -9.32483  ]
 2 : [  1.40226  ]
 3 : [ -7.32298  ]
 4 : [ -0.218334 ]
 5 : [  2.36082  ]
 6 : [ -0.995967 ]
 7 : [ -4.08345  ]
 8 : [ -6.88382  ]
 9 : [ -9.71897  ]
10 : [ -6.76925  ]
11 : [ -2.46257  ]
Points NOT verifying the Event outside the vicinity of the design points =
Model values =