otbenchmark
Note
Go to the end to download the full example code.
Check reliability problems reference probabilities¶
In this example, we check that the reference probabilities in the reliability problems are consistent with confidence bounds from Monte-Carlo simulations. These 95% confidence bounds are stored in ‘reliability_compute_reference_proba.csv’ and required approximately than :math`10^9` function evaluations for each problem.
We consider two different metrics:
we check if the reference probability is within the 95% confidence bounds,
we compute the number of significant digits by comparing the Monte-Carlo estimator and the reference value.
The number of significant digits may be as high as 17 when all decimal digits are correct. However, the reference probabilities are only known up to 3 digits for most problems. In order to keep some safeguard, we will be happy with 2 correct digits.
These metrics may fail.
On average, a fraction equal to 5% of the estimates are outside the confidence bounds.
The Monte-Carlo estimator may not be accurate enough, e.g. if the probability is very close to zero.
import otbenchmark as otb
import pandas as pd
df = pd.read_csv("reliability_compute_reference_proba.csv")
df
data = df.values
pf_reference = data[:, 1]
pmin = data[:, 3]
pmax = data[:, 4]
benchmarkProblemList = otb.ReliabilityBenchmarkProblemList()
numberOfProblems = len(benchmarkProblemList)
numberOfProblems
26
digitsMinimum = 2
categoryA = []
categoryB = []
categoryC = []
categoryD = []
for i in range(numberOfProblems):
problem = benchmarkProblemList[i]
name = problem.getName()
pf = problem.getProbability()
event = problem.getEvent()
antecedent = event.getAntecedent()
distribution = antecedent.getDistribution()
dimension = distribution.getDimension()
if pf > pmin[i] and pf < pmax[i]:
tagBounds = "In"
else:
tagBounds = "Out"
digits = otb.ComputeLogRelativeError(pf_reference[i], pf)
if tagBounds == "In" and digits >= digitsMinimum:
categoryA.append(name)
elif tagBounds == "Out" and digits >= digitsMinimum:
categoryB.append(name)
elif tagBounds == "In" and digits < digitsMinimum:
categoryC.append(name)
else:
categoryD.append(name)
print(
"#%d, %-10s, pf=%.2e, ref=%.2e, C.I.=[%.2e,%.2e], digits=%d : %s"
% (i, name[0:10], pf, pf_reference[i], pmin[i], pmax[i], digits, tagBounds)
)
#0, RP8 , pf=7.90e-04, ref=7.91e-04, C.I.=[7.87e-04,7.94e-04], digits=2 : In
#1, RP14 , pf=7.73e-04, ref=7.71e-04, C.I.=[7.69e-04,7.73e-04], digits=2 : In
#2, RP22 , pf=4.21e-03, ref=4.21e-03, C.I.=[4.20e-03,4.21e-03], digits=4 : In
#3, RP24 , pf=2.86e-03, ref=2.86e-03, C.I.=[2.86e-03,2.86e-03], digits=3 : In
#4, RP25 , pf=4.15e-05, ref=4.18e-05, C.I.=[4.14e-05,4.21e-05], digits=2 : In
#5, RP28 , pf=1.45e-07, ref=1.32e-07, C.I.=[1.15e-07,1.48e-07], digits=0 : In
#6, RP31 , pf=3.23e-03, ref=3.23e-03, C.I.=[3.22e-03,3.23e-03], digits=3 : In
#7, RP33 , pf=2.57e-03, ref=2.57e-03, C.I.=[2.57e-03,2.58e-03], digits=2 : Out
#8, RP35 , pf=3.48e-03, ref=3.48e-03, C.I.=[3.48e-03,3.48e-03], digits=5 : In
#9, RP38 , pf=8.10e-03, ref=8.06e-03, C.I.=[8.05e-03,8.07e-03], digits=2 : Out
#10, RP53 , pf=3.13e-02, ref=3.13e-02, C.I.=[3.13e-02,3.13e-02], digits=3 : Out
#11, RP55 , pf=5.60e-01, ref=5.60e-01, C.I.=[5.60e-01,5.60e-01], digits=4 : In
#12, RP54 , pf=9.98e-04, ref=9.93e-04, C.I.=[9.88e-04,9.97e-04], digits=2 : Out
#13, RP57 , pf=2.84e-02, ref=2.82e-02, C.I.=[2.82e-02,2.82e-02], digits=2 : Out
#14, RP75 , pf=9.82e-03, ref=9.82e-03, C.I.=[9.81e-03,9.82e-03], digits=4 : In
#15, RP89 , pf=5.43e-03, ref=5.47e-03, C.I.=[5.47e-03,5.47e-03], digits=2 : Out
#16, RP107 , pf=2.92e-07, ref=2.75e-07, C.I.=[2.30e-07,3.20e-07], digits=1 : In
#17, RP110 , pf=3.19e-05, ref=3.18e-05, C.I.=[3.15e-05,3.21e-05], digits=2 : In
#18, RP111 , pf=7.65e-07, ref=7.85e-07, C.I.=[7.41e-07,8.29e-07], digits=1 : In
#19, RP63 , pf=3.79e-04, ref=3.77e-04, C.I.=[3.72e-04,3.82e-04], digits=2 : In
#20, RP91 , pf=6.97e-04, ref=7.00e-04, C.I.=[6.98e-04,7.02e-04], digits=2 : Out
#21, RP60 , pf=4.56e-02, ref=4.48e-02, C.I.=[4.48e-02,4.49e-02], digits=1 : Out
#22, RP77 , pf=2.87e-07, ref=2.71e-07, C.I.=[2.43e-07,2.99e-07], digits=1 : In
#23, Four-branc, pf=2.22e-03, ref=2.23e-03, C.I.=[2.22e-03,2.23e-03], digits=2 : In
#24, R-S , pf=7.86e-02, ref=7.86e-02, C.I.=[7.86e-02,7.87e-02], digits=4 : In
#25, Axial stre, pf=2.92e-02, ref=2.92e-02, C.I.=[2.92e-02,2.92e-02], digits=4 : In
There are four different cases.
Category A: all good. For some problems, both metrics are correct in the sense that the reference probability is within the bounds and the number of significant digits is larger than 2. The RP24, RP55, RP110, RP63, R-S, Axial stressed beam problems fall in that category.
Category B: correct digits, not in bounds. We see that the RP8 problem has a reference probability outside of the 95% confidence bounds, but has 2 significant digits.
Category C: insufficient digits, in bounds. The difficult RP28 problem fall in that category.
Category D: insufficient digits, not in bounds. These are suspicious problems.
print(categoryA)
['RP8', 'RP14', 'RP22', 'RP24', 'RP25', 'RP31', 'RP35', 'RP55', 'RP75', 'RP110', 'RP63', 'Four-branch serial system', 'R-S', 'Axial stressed beam']
print(categoryB)
['RP33', 'RP38', 'RP53', 'RP54', 'RP57', 'RP89', 'RP91']
print(categoryC)
['RP28', 'RP107', 'RP111', 'RP77']
print(categoryD)
['RP60']
The number of suspicious problems seems very large. However, we notice that all these cases are so that the reference probability is close, in absolute value, to the Monte-Carlo estimator.
Total running time of the script: (0 minutes 0.498 seconds)