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 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
Unnamed: 0 PF N. function calls PMin PMax C.O.V. Digits Time (s)
0 RP8 7.908179e-04 2.413400e+08 7.872713e-04 7.943646e-04 0.002288 1.640503 300.011341
1 RP14 7.708905e-04 7.441900e+08 7.688964e-04 7.728846e-04 0.001320 1.879484 300.004585
2 RP22 4.207357e-03 1.495610e+09 4.204076e-03 4.210637e-03 0.000398 2.400308 300.000947
3 RP24 2.860848e-03 1.644700e+09 2.858266e-03 2.863429e-03 0.000460 2.336891 300.001723
4 RP25 4.175883e-05 1.567860e+09 4.143895e-05 4.207871e-05 0.003908 1.408015 300.015279
5 RP28 1.315725e-07 1.839290e+09 1.149947e-07 1.481503e-07 0.064286 0.191886 300.000353
6 RP31 3.227556e-03 1.787260e+09 3.224926e-03 3.230186e-03 0.000416 2.381211 300.000409
7 RP33 2.574817e-03 1.430700e+09 2.572190e-03 2.577443e-03 0.000520 2.283686 300.000557
8 RP35 3.478964e-03 1.415140e+09 3.475896e-03 3.482032e-03 0.000450 2.346860 300.002611
9 RP38 8.059349e-03 7.765700e+08 8.053061e-03 8.065638e-03 0.000398 2.399976 300.003572
10 RP53 3.131966e-02 1.420390e+09 3.131060e-02 3.132872e-02 0.000148 2.831000 300.002226
11 RP55 5.600269e-01 1.523630e+09 5.600020e-01 5.600519e-01 0.000023 3.643809 300.000719
12 RP54 9.927480e-04 1.715400e+08 9.880351e-04 9.974610e-04 0.002422 1.615796 300.018069
13 RP57 2.822772e-02 1.249510e+09 2.821854e-02 2.823691e-02 0.000166 2.779904 300.000954
14 RP75 9.818417e-03 1.597510e+09 9.813582e-03 9.823253e-03 0.000251 2.599863 300.000640
15 RP89 5.469847e-03 1.366180e+09 5.465935e-03 5.473758e-03 0.000365 2.437911 300.000577
16 RP107 2.754926e-07 5.227000e+08 2.304941e-07 3.204912e-07 0.083337 0.079160 300.004477
17 RP110 3.183607e-05 1.344010e+09 3.153441e-05 3.213774e-05 0.004835 1.315646 300.001033
18 RP111 7.851043e-07 1.552660e+09 7.410290e-07 8.291796e-07 0.028643 0.542980 300.000674
19 RP63 3.772015e-04 5.360000e+07 3.720028e-04 3.824002e-04 0.007032 1.152929 300.044605
20 RP91 6.998219e-04 8.196500e+08 6.980114e-04 7.016324e-04 0.001320 1.879438 300.002981
21 RP60 4.483579e-02 3.014400e+08 4.481243e-02 4.485916e-02 0.000266 2.575352 300.002735
22 RP77 2.711905e-07 1.342230e+09 2.433297e-07 2.990513e-07 0.052417 0.280529 300.000783
23 Four-branch serial system 2.225032e-03 1.351070e+09 2.222519e-03 2.227545e-03 0.000576 2.239469 300.000255
24 R-S 7.864349e-02 1.839670e+09 7.863119e-02 7.865579e-02 0.000080 3.097966 300.000382
25 Axial stressed beam 2.919903e-02 1.391830e+09 2.919019e-02 2.920788e-02 0.000155 2.810891 300.001095


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.430 seconds)