Compute reference probabilities with Monte-Carlo

In this example, we perform a reliability benchmark based on a large Monte-Carlo sample. In order to limit the elapsed time, we consider a limited elapsed time for each problem. In order to get the best possible accuracy within this time limit, we set the coefficient of variation to zero.

import otbenchmark as otb
import pandas as pd
import numpy as np
from tqdm import tqdm
import time

problemslist = otb.ReliabilityBenchmarkProblemList()
numberOfProblems = len(problemslist)
numberOfProblems

26

coefficientOfVariation = 0.0
maximumOuterSampling = 10 ** 4  # 10 ** 6 for real
blockSize = 10 ** 0  # 10 ** 4 for real simulations
blockSize

1

confidenceLevel = 0.95
maximumDurationSeconds = 5 * 60.0

totalDurationMinutes = numberOfProblems * maximumDurationSeconds / 60.0
totalDurationMinutes

130.0

model_names = [problemslist[i].getName() for i in range(numberOfProblems)]
metrics = ["PF", "N. function calls", "PMin", "PMax", "C.O.V.", "Digits", "Time (s)"]
resultArray = np.zeros((numberOfProblems, len(metrics)))
for i in tqdm(range(numberOfProblems)):
    startTime = time.time()
    problem = problemslist[i]
    name = problem.getName()
    event = problem.getEvent()
    g = event.getFunction()
    factory = otb.ProbabilitySimulationAlgorithmFactory()
    algo = factory.buildMonteCarlo(problem)
    algo.setMaximumOuterSampling(maximumOuterSampling)
    algo.setBlockSize(blockSize)
    algo.setMaximumCoefficientOfVariation(coefficientOfVariation)
    algo.setMaximumTimeDuration(maximumDurationSeconds)
    initialNumberOfCall = g.getEvaluationCallsNumber()
    algo.run()
    result = algo.getResult()
    numberOfFunctionEvaluations = g.getEvaluationCallsNumber() - initialNumberOfCall
    computedProbability = result.getProbabilityEstimate()
    confidenceLength = result.getConfidenceLength(confidenceLevel)
    pmin = computedProbability - 0.5 * confidenceLength
    pmax = computedProbability + 0.5 * confidenceLength
    cov = result.getCoefficientOfVariation()
    if cov > 0.0:
        expectedDigits = -np.log10(cov) - 1.0
    else:
        expectedDigits = 0.0
    stopTime = time.time()
    elapsedTime = stopTime - startTime
    resultArray[i][0] = computedProbability
    resultArray[i][1] = numberOfFunctionEvaluations
    resultArray[i][2] = pmin
    resultArray[i][3] = pmax
    resultArray[i][4] = cov
    resultArray[i][5] = expectedDigits
    resultArray[i][6] = elapsedTime

  0%|          | 0/26 [00:00<?, ?it/s]
  4%|▍         | 1/26 [00:00<00:08,  3.12it/s]
 12%|█▏        | 3/26 [00:00<00:02,  8.04it/s]
 31%|███       | 8/26 [00:00<00:00, 19.07it/s]
 46%|████▌     | 12/26 [00:00<00:00, 23.64it/s]
 58%|█████▊    | 15/26 [00:00<00:00, 23.41it/s]
 69%|██████▉   | 18/26 [00:00<00:00, 24.33it/s]
 81%|████████  | 21/26 [00:01<00:00, 12.07it/s]
 88%|████████▊ | 23/26 [00:01<00:00, 10.14it/s]
100%|██████████| 26/26 [00:01<00:00, 14.27it/s]

df = pd.DataFrame(resultArray, index=model_names, columns=metrics)
# df.to_csv("reliability_compute_reference_proba.csv")
df

	PF	N. function calls	PMin	PMax	C.O.V.	Digits	Time (s)
RP8	0.0010	10000.0	0.000381	0.001619	0.316070	-0.499783	0.320407
RP14	0.0007	10000.0	0.000182	0.001218	0.377832	-0.577299	0.087304
RP22	0.0044	10000.0	0.003103	0.005697	0.150424	-0.177316	0.023932
RP24	0.0028	10000.0	0.001764	0.003836	0.188717	-0.275812	0.023893
RP25	0.0001	10000.0	-0.000096	0.000296	0.999950	-0.999978	0.024736
RP28	0.0000	10000.0	0.000000	0.000000	-1.000000	0.000000	0.025470
RP31	0.0034	10000.0	0.002259	0.004541	0.171207	-0.233521	0.023821
RP33	0.0029	10000.0	0.001846	0.003954	0.185426	-0.268170	0.028183
RP35	0.0037	10000.0	0.002510	0.004890	0.164095	-0.215094	0.024314
RP38	0.0077	10000.0	0.005987	0.009413	0.113521	-0.055076	0.046152
RP53	0.0342	10000.0	0.030638	0.037762	0.053141	0.274569	0.024471
RP55	0.5565	10000.0	0.546763	0.566237	0.008927	1.049286	0.025169
RP54	0.0009	10000.0	0.000312	0.001488	0.333183	-0.522683	0.082615
RP57	0.0262	10000.0	0.023069	0.029331	0.060966	0.214916	0.024346
RP75	0.0098	10000.0	0.007869	0.011731	0.100519	-0.002248	0.023711
RP89	0.0051	10000.0	0.003704	0.006496	0.139670	-0.145105	0.024375
RP107	0.0000	10000.0	0.000000	0.000000	-1.000000	0.000000	0.059876
RP110	0.0000	10000.0	0.000000	0.000000	-1.000000	0.000000	0.028418
RP111	0.0000	10000.0	0.000000	0.000000	-1.000000	0.000000	0.036227
RP63	0.0003	10000.0	-0.000039	0.000639	0.577264	-0.761374	0.448361
RP91	0.0006	10000.0	0.000120	0.001080	0.408126	-0.610794	0.038954
RP60	0.0426	10000.0	0.038642	0.046558	0.047407	0.324158	0.272267
RP77	0.0000	10000.0	0.000000	0.000000	-1.000000	0.000000	0.030778
Four-branch serial system	0.0019	10000.0	0.001046	0.002754	0.229198	-0.360210	0.024731
R-S	0.0791	10000.0	0.073810	0.084390	0.034121	0.466982	0.023885
Axial stressed beam	0.0277	10000.0	0.024483	0.030917	0.059246	0.227340	0.024011

The problems with higher failture probabilities are obviously solved with more accuracy with the Monte-Carlo method. For example, the RP55 problem which has the highest probability equal to 0.560 has more than 3 significant digits. On the opposite side, the problems with probabilities close to zero are much more difficult to solve. The RP28 with a probability close to has no significant digit.

These previous results are consistent with the distribution of the Monte-Carlo estimator. The properties of the binomial distribution imply that its variance is:

where is the sample size and is the failure probability. The coefficient of variation is:

Since we do not know the exact value of , we use is approximation instead (this turns rigorous equations into approximate ones, but does not change the outcome). This implies:

When , we have which implies:

Inverting the previous equation, we get the sample size given the coefficient of variation:

This leads to the rule of thumb that, in order to estimate the probability , where is an integer, we need a sample size equal to:

For example, estimating the probability of the RP28 problem with just one single digit leads to a sample size equal to , since the exact .