Convergence of Monte-Carlo to estimate the probability in a reliability problem¶

The goal of this document is to present the convergence of the Monte-Carlo algorithm to the exact probability when the sample size increases. This convergence is expressed in terms of absolute error. We show that the rate of convergence is $O(\sqrt{n})$ , where $n$ is the sample size.

import openturns as ot
import openturns.viewer as otv
import numpy as np
import otbenchmark as otb
import time

problem = otb.RminusSReliability()
print(problem)


def ComputeProbabilityFromMonteCarlo(
    problem, coefficientOfVariation=0.1, maximumOuterSampling=1000, blockSize=2
):
    event = problem.getEvent()
    g = event.getFunction()
    # Create the Monte-Carlo algorithm
    algoProb = ot.ProbabilitySimulationAlgorithm(event)
    algoProb.setMaximumOuterSampling(maximumOuterSampling)
    algoProb.setBlockSize(blockSize)
    algoProb.setMaximumCoefficientOfVariation(coefficientOfVariation)
    initialNumberOfFunctionEvaluations = g.getEvaluationCallsNumber()
    algoProb.run()
    # Get the results
    resultAlgo = algoProb.getResult()
    numberOfFunctionEvaluations = (
        g.getEvaluationCallsNumber() - initialNumberOfFunctionEvaluations
    )
    pf = resultAlgo.getProbabilityEstimate()
    level = 0.95
    c95 = resultAlgo.getConfidenceLength(level)
    pmin = pf - 0.5 * c95
    pmax = pf + 0.5 * c95
    print(
        "Number of function calls = %d" % (numberOfFunctionEvaluations),
        ", Pf = %.4f" % (pf),
        ", %.1f %% confidence interval :[%.4f,%.4f] " % (level * 100, pmin, pmax),
    )
    absoluteError = abs(pf - problem.getProbability())
    result = {
        "numberOfFunctionEvaluations": numberOfFunctionEvaluations,
        "pf": pf,
        "pmin": pmin,
        "pmax": pmax,
        "absoluteError": absoluteError,
    }
    return result

name = R-S
event = class=ThresholdEventImplementation antecedent=class=CompositeRandomVector function=class=Function name=Unnamed implementation=class=FunctionImplementation name=Unnamed description=[R,S,y0] evaluationImplementation=class=SymbolicEvaluation name=Unnamed inputVariablesNames=[R,S] outputVariablesNames=[y0] formulas=[R - S] gradientImplementation=class=SymbolicGradient name=Unnamed evaluation=class=SymbolicEvaluation name=Unnamed inputVariablesNames=[R,S] outputVariablesNames=[y0] formulas=[R - S] hessianImplementation=class=SymbolicHessian name=Unnamed evaluation=class=SymbolicEvaluation name=Unnamed inputVariablesNames=[R,S] outputVariablesNames=[y0] formulas=[R - S] antecedent=class=UsualRandomVector distribution=class=JointDistribution name=JointDistribution dimension=2 copula=class=IndependentCopula name=IndependentCopula dimension=2 marginal[0]=class=Normal name=Normal dimension=1 mean=class=Point name=Unnamed dimension=1 values=[4] sigma=class=Point name=Unnamed dimension=1 values=[1] correlationMatrix=class=CorrelationMatrix dimension=1 implementation=class=MatrixImplementation name=Unnamed rows=1 columns=1 values=[1] marginal[1]=class=Normal name=Normal dimension=1 mean=class=Point name=Unnamed dimension=1 values=[2] sigma=class=Point name=Unnamed dimension=1 values=[1] correlationMatrix=class=CorrelationMatrix dimension=1 implementation=class=MatrixImplementation name=Unnamed rows=1 columns=1 values=[1] operator=class=Less name=Unnamed threshold=0
probability = 0.07864960352514257

result = ComputeProbabilityFromMonteCarlo(problem)

Number of function calls = 1978 , Pf = 0.0713 , 95.0 % confidence interval :[0.0573,0.0852]

numberOfPoints = 15  # Number of atomic experiments
numberOfRepetitions = 10  # Number of repetitions of each experiment
sampleSizeAbsoluteErrorTable = ot.Sample(numberOfPoints * numberOfRepetitions, 2)
sampleSizeAbsoluteErrorTable.setDescription(["Sample size", "Absolute error"])

cov = 0.0
startTime = time.time()

maximumOuterSampling = 1
index = 0
for i in range(numberOfPoints):
    maximumOuterSampling *= 2
    for j in range(numberOfRepetitions):
        result = ComputeProbabilityFromMonteCarlo(
            problem,
            coefficientOfVariation=cov,
            maximumOuterSampling=maximumOuterSampling,
        )
        sampleSizeAbsoluteErrorTable[index, 0] = result["numberOfFunctionEvaluations"]
        sampleSizeAbsoluteErrorTable[index, 1] = result["absoluteError"]
        index += 1

Number of function calls = 4 , Pf = 0.0000 , 95.0 % confidence interval :[0.0000,0.0000]
Number of function calls = 4 , Pf = 0.0000 , 95.0 % confidence interval :[0.0000,0.0000]
Number of function calls = 4 , Pf = 0.2500 , 95.0 % confidence interval :[-0.2978,0.7978]
Number of function calls = 4 , Pf = 0.5000 , 95.0 % confidence interval :[0.0100,0.9900]
Number of function calls = 4 , Pf = 0.0000 , 95.0 % confidence interval :[0.0000,0.0000]
Number of function calls = 4 , Pf = 0.0000 , 95.0 % confidence interval :[0.0000,0.0000]
Number of function calls = 4 , Pf = 0.0000 , 95.0 % confidence interval :[0.0000,0.0000]
Number of function calls = 4 , Pf = 0.0000 , 95.0 % confidence interval :[0.0000,0.0000]
Number of function calls = 4 , Pf = 0.0000 , 95.0 % confidence interval :[0.0000,0.0000]
Number of function calls = 4 , Pf = 0.0000 , 95.0 % confidence interval :[0.0000,0.0000]
Number of function calls = 8 , Pf = 0.0000 , 95.0 % confidence interval :[0.0000,0.0000]
Number of function calls = 8 , Pf = 0.1250 , 95.0 % confidence interval :[-0.1623,0.4123]
Number of function calls = 8 , Pf = 0.0000 , 95.0 % confidence interval :[0.0000,0.0000]
Number of function calls = 8 , Pf = 0.1250 , 95.0 % confidence interval :[-0.1623,0.4123]
Number of function calls = 8 , Pf = 0.3750 , 95.0 % confidence interval :[-0.0026,0.7526]
Number of function calls = 8 , Pf = 0.2500 , 95.0 % confidence interval :[-0.0501,0.5501]
Number of function calls = 8 , Pf = 0.0000 , 95.0 % confidence interval :[0.0000,0.0000]
Number of function calls = 8 , Pf = 0.2500 , 95.0 % confidence interval :[-0.1374,0.6374]
Number of function calls = 8 , Pf = 0.0000 , 95.0 % confidence interval :[0.0000,0.0000]
Number of function calls = 8 , Pf = 0.1250 , 95.0 % confidence interval :[-0.1623,0.4123]
Number of function calls = 16 , Pf = 0.1250 , 95.0 % confidence interval :[-0.0781,0.3281]
Number of function calls = 16 , Pf = 0.0625 , 95.0 % confidence interval :[-0.0844,0.2094]
Number of function calls = 16 , Pf = 0.0625 , 95.0 % confidence interval :[-0.0844,0.2094]
Number of function calls = 16 , Pf = 0.0625 , 95.0 % confidence interval :[-0.0844,0.2094]
Number of function calls = 16 , Pf = 0.0625 , 95.0 % confidence interval :[-0.0844,0.2094]
Number of function calls = 16 , Pf = 0.1250 , 95.0 % confidence interval :[-0.0370,0.2870]
Number of function calls = 16 , Pf = 0.0625 , 95.0 % confidence interval :[-0.0844,0.2094]
Number of function calls = 16 , Pf = 0.0625 , 95.0 % confidence interval :[-0.0844,0.2094]
Number of function calls = 16 , Pf = 0.1250 , 95.0 % confidence interval :[-0.0781,0.3281]
Number of function calls = 16 , Pf = 0.0000 , 95.0 % confidence interval :[0.0000,0.0000]
Number of function calls = 32 , Pf = 0.0312 , 95.0 % confidence interval :[-0.0430,0.1055]
Number of function calls = 32 , Pf = 0.0938 , 95.0 % confidence interval :[-0.0321,0.2196]
Number of function calls = 32 , Pf = 0.0625 , 95.0 % confidence interval :[-0.0414,0.1664]
Number of function calls = 32 , Pf = 0.0938 , 95.0 % confidence interval :[-0.0321,0.2196]
Number of function calls = 32 , Pf = 0.0000 , 95.0 % confidence interval :[0.0000,0.0000]
Number of function calls = 32 , Pf = 0.1875 , 95.0 % confidence interval :[0.0390,0.3360]
Number of function calls = 32 , Pf = 0.0938 , 95.0 % confidence interval :[-0.0321,0.2196]
Number of function calls = 32 , Pf = 0.1562 , 95.0 % confidence interval :[-0.0025,0.3150]
Number of function calls = 32 , Pf = 0.0625 , 95.0 % confidence interval :[-0.0414,0.1664]
Number of function calls = 32 , Pf = 0.1250 , 95.0 % confidence interval :[-0.0186,0.2686]
Number of function calls = 64 , Pf = 0.0781 , 95.0 % confidence interval :[-0.0035,0.1598]
Number of function calls = 64 , Pf = 0.0938 , 95.0 % confidence interval :[0.0048,0.1827]
Number of function calls = 64 , Pf = 0.0469 , 95.0 % confidence interval :[-0.0171,0.1108]
Number of function calls = 64 , Pf = 0.1562 , 95.0 % confidence interval :[0.0440,0.2685]
Number of function calls = 64 , Pf = 0.0938 , 95.0 % confidence interval :[0.0048,0.1827]
Number of function calls = 64 , Pf = 0.0469 , 95.0 % confidence interval :[-0.0171,0.1108]
Number of function calls = 64 , Pf = 0.1719 , 95.0 % confidence interval :[0.0589,0.2849]
Number of function calls = 64 , Pf = 0.0781 , 95.0 % confidence interval :[0.0024,0.1538]
Number of function calls = 64 , Pf = 0.0938 , 95.0 % confidence interval :[0.0102,0.1773]
Number of function calls = 64 , Pf = 0.0469 , 95.0 % confidence interval :[-0.0171,0.1108]
Number of function calls = 128 , Pf = 0.1094 , 95.0 % confidence interval :[0.0418,0.1769]
Number of function calls = 128 , Pf = 0.0312 , 95.0 % confidence interval :[-0.0059,0.0684]
Number of function calls = 128 , Pf = 0.1250 , 95.0 % confidence interval :[0.0548,0.1952]
Number of function calls = 128 , Pf = 0.0312 , 95.0 % confidence interval :[-0.0059,0.0684]
Number of function calls = 128 , Pf = 0.0469 , 95.0 % confidence interval :[0.0017,0.0921]
Number of function calls = 128 , Pf = 0.0391 , 95.0 % confidence interval :[-0.0023,0.0804]
Number of function calls = 128 , Pf = 0.0859 , 95.0 % confidence interval :[0.0275,0.1444]
Number of function calls = 128 , Pf = 0.1172 , 95.0 % confidence interval :[0.0492,0.1852]
Number of function calls = 128 , Pf = 0.0703 , 95.0 % confidence interval :[0.0154,0.1252]
Number of function calls = 128 , Pf = 0.1016 , 95.0 % confidence interval :[0.0381,0.1650]
Number of function calls = 256 , Pf = 0.0820 , 95.0 % confidence interval :[0.0417,0.1224]
Number of function calls = 256 , Pf = 0.0859 , 95.0 % confidence interval :[0.0432,0.1286]
Number of function calls = 256 , Pf = 0.0781 , 95.0 % confidence interval :[0.0373,0.1190]
Number of function calls = 256 , Pf = 0.0625 , 95.0 % confidence interval :[0.0266,0.0984]
Number of function calls = 256 , Pf = 0.0781 , 95.0 % confidence interval :[0.0380,0.1182]
Number of function calls = 256 , Pf = 0.0352 , 95.0 % confidence interval :[0.0074,0.0630]
Number of function calls = 256 , Pf = 0.0664 , 95.0 % confidence interval :[0.0302,0.1026]
Number of function calls = 256 , Pf = 0.0547 , 95.0 % confidence interval :[0.0202,0.0891]
Number of function calls = 256 , Pf = 0.0625 , 95.0 % confidence interval :[0.0258,0.0992]
Number of function calls = 256 , Pf = 0.0703 , 95.0 % confidence interval :[0.0315,0.1092]
Number of function calls = 512 , Pf = 0.0723 , 95.0 % confidence interval :[0.0450,0.0996]
Number of function calls = 512 , Pf = 0.0605 , 95.0 % confidence interval :[0.0358,0.0852]
Number of function calls = 512 , Pf = 0.0996 , 95.0 % confidence interval :[0.0686,0.1306]
Number of function calls = 512 , Pf = 0.0938 , 95.0 % confidence interval :[0.0628,0.1247]
Number of function calls = 512 , Pf = 0.0742 , 95.0 % confidence interval :[0.0471,0.1013]
Number of function calls = 512 , Pf = 0.0723 , 95.0 % confidence interval :[0.0450,0.0996]
Number of function calls = 512 , Pf = 0.0684 , 95.0 % confidence interval :[0.0415,0.0952]
Number of function calls = 512 , Pf = 0.0859 , 95.0 % confidence interval :[0.0565,0.1154]
Number of function calls = 512 , Pf = 0.0938 , 95.0 % confidence interval :[0.0628,0.1247]
Number of function calls = 512 , Pf = 0.0566 , 95.0 % confidence interval :[0.0322,0.0811]
Number of function calls = 1024 , Pf = 0.0713 , 95.0 % confidence interval :[0.0519,0.0906]
Number of function calls = 1024 , Pf = 0.0791 , 95.0 % confidence interval :[0.0588,0.0994]
Number of function calls = 1024 , Pf = 0.0850 , 95.0 % confidence interval :[0.0638,0.1061]
Number of function calls = 1024 , Pf = 0.0850 , 95.0 % confidence interval :[0.0640,0.1059]
Number of function calls = 1024 , Pf = 0.0830 , 95.0 % confidence interval :[0.0625,0.1035]
Number of function calls = 1024 , Pf = 0.0869 , 95.0 % confidence interval :[0.0658,0.1080]
Number of function calls = 1024 , Pf = 0.0869 , 95.0 % confidence interval :[0.0657,0.1081]
Number of function calls = 1024 , Pf = 0.0771 , 95.0 % confidence interval :[0.0570,0.0973]
Number of function calls = 1024 , Pf = 0.0830 , 95.0 % confidence interval :[0.0622,0.1038]
Number of function calls = 1024 , Pf = 0.0742 , 95.0 % confidence interval :[0.0546,0.0939]
Number of function calls = 2048 , Pf = 0.0840 , 95.0 % confidence interval :[0.0693,0.0987]
Number of function calls = 2048 , Pf = 0.0864 , 95.0 % confidence interval :[0.0716,0.1013]
Number of function calls = 2048 , Pf = 0.0894 , 95.0 % confidence interval :[0.0743,0.1044]
Number of function calls = 2048 , Pf = 0.0864 , 95.0 % confidence interval :[0.0715,0.1014]
Number of function calls = 2048 , Pf = 0.0757 , 95.0 % confidence interval :[0.0618,0.0896]
Number of function calls = 2048 , Pf = 0.0781 , 95.0 % confidence interval :[0.0638,0.0925]
Number of function calls = 2048 , Pf = 0.0786 , 95.0 % confidence interval :[0.0643,0.0930]
Number of function calls = 2048 , Pf = 0.0747 , 95.0 % confidence interval :[0.0609,0.0885]
Number of function calls = 2048 , Pf = 0.0674 , 95.0 % confidence interval :[0.0542,0.0806]
Number of function calls = 2048 , Pf = 0.0796 , 95.0 % confidence interval :[0.0653,0.0938]
Number of function calls = 4096 , Pf = 0.0735 , 95.0 % confidence interval :[0.0637,0.0833]
Number of function calls = 4096 , Pf = 0.0774 , 95.0 % confidence interval :[0.0674,0.0874]
Number of function calls = 4096 , Pf = 0.0825 , 95.0 % confidence interval :[0.0722,0.0929]
Number of function calls = 4096 , Pf = 0.0764 , 95.0 % confidence interval :[0.0665,0.0863]
Number of function calls = 4096 , Pf = 0.0815 , 95.0 % confidence interval :[0.0713,0.0918]
Number of function calls = 4096 , Pf = 0.0791 , 95.0 % confidence interval :[0.0689,0.0893]
Number of function calls = 4096 , Pf = 0.0784 , 95.0 % confidence interval :[0.0682,0.0885]
Number of function calls = 4096 , Pf = 0.0740 , 95.0 % confidence interval :[0.0641,0.0838]
Number of function calls = 4096 , Pf = 0.0747 , 95.0 % confidence interval :[0.0649,0.0845]
Number of function calls = 4096 , Pf = 0.0815 , 95.0 % confidence interval :[0.0712,0.0918]
Number of function calls = 8192 , Pf = 0.0782 , 95.0 % confidence interval :[0.0711,0.0854]
Number of function calls = 8192 , Pf = 0.0781 , 95.0 % confidence interval :[0.0710,0.0852]
Number of function calls = 8192 , Pf = 0.0770 , 95.0 % confidence interval :[0.0700,0.0841]
Number of function calls = 8192 , Pf = 0.0768 , 95.0 % confidence interval :[0.0697,0.0838]
Number of function calls = 8192 , Pf = 0.0745 , 95.0 % confidence interval :[0.0675,0.0815]
Number of function calls = 8192 , Pf = 0.0780 , 95.0 % confidence interval :[0.0709,0.0851]
Number of function calls = 8192 , Pf = 0.0797 , 95.0 % confidence interval :[0.0725,0.0869]
Number of function calls = 8192 , Pf = 0.0797 , 95.0 % confidence interval :[0.0725,0.0869]
Number of function calls = 8192 , Pf = 0.0760 , 95.0 % confidence interval :[0.0690,0.0831]
Number of function calls = 8192 , Pf = 0.0818 , 95.0 % confidence interval :[0.0745,0.0891]
Number of function calls = 16384 , Pf = 0.0825 , 95.0 % confidence interval :[0.0773,0.0876]
Number of function calls = 16384 , Pf = 0.0812 , 95.0 % confidence interval :[0.0761,0.0863]
Number of function calls = 16384 , Pf = 0.0760 , 95.0 % confidence interval :[0.0710,0.0810]
Number of function calls = 16384 , Pf = 0.0790 , 95.0 % confidence interval :[0.0740,0.0841]
Number of function calls = 16384 , Pf = 0.0793 , 95.0 % confidence interval :[0.0743,0.0844]
Number of function calls = 16384 , Pf = 0.0806 , 95.0 % confidence interval :[0.0755,0.0857]
Number of function calls = 16384 , Pf = 0.0785 , 95.0 % confidence interval :[0.0735,0.0835]
Number of function calls = 16384 , Pf = 0.0779 , 95.0 % confidence interval :[0.0729,0.0830]
Number of function calls = 16384 , Pf = 0.0786 , 95.0 % confidence interval :[0.0736,0.0837]
Number of function calls = 16384 , Pf = 0.0750 , 95.0 % confidence interval :[0.0701,0.0799]
Number of function calls = 32768 , Pf = 0.0772 , 95.0 % confidence interval :[0.0737,0.0808]
Number of function calls = 32768 , Pf = 0.0766 , 95.0 % confidence interval :[0.0731,0.0801]
Number of function calls = 32768 , Pf = 0.0796 , 95.0 % confidence interval :[0.0760,0.0832]
Number of function calls = 32768 , Pf = 0.0802 , 95.0 % confidence interval :[0.0766,0.0838]
Number of function calls = 32768 , Pf = 0.0796 , 95.0 % confidence interval :[0.0760,0.0831]
Number of function calls = 32768 , Pf = 0.0772 , 95.0 % confidence interval :[0.0737,0.0807]
Number of function calls = 32768 , Pf = 0.0786 , 95.0 % confidence interval :[0.0750,0.0822]
Number of function calls = 32768 , Pf = 0.0778 , 95.0 % confidence interval :[0.0743,0.0814]
Number of function calls = 32768 , Pf = 0.0810 , 95.0 % confidence interval :[0.0774,0.0846]
Number of function calls = 32768 , Pf = 0.0808 , 95.0 % confidence interval :[0.0772,0.0844]
Number of function calls = 65536 , Pf = 0.0790 , 95.0 % confidence interval :[0.0765,0.0815]
Number of function calls = 65536 , Pf = 0.0774 , 95.0 % confidence interval :[0.0749,0.0799]
Number of function calls = 65536 , Pf = 0.0786 , 95.0 % confidence interval :[0.0761,0.0812]
Number of function calls = 65536 , Pf = 0.0781 , 95.0 % confidence interval :[0.0756,0.0807]
Number of function calls = 65536 , Pf = 0.0793 , 95.0 % confidence interval :[0.0767,0.0818]
Number of function calls = 65536 , Pf = 0.0793 , 95.0 % confidence interval :[0.0768,0.0818]
Number of function calls = 65536 , Pf = 0.0806 , 95.0 % confidence interval :[0.0781,0.0832]
Number of function calls = 65536 , Pf = 0.0784 , 95.0 % confidence interval :[0.0759,0.0810]
Number of function calls = 65536 , Pf = 0.0783 , 95.0 % confidence interval :[0.0758,0.0808]
Number of function calls = 65536 , Pf = 0.0789 , 95.0 % confidence interval :[0.0763,0.0814]

elapsedTime = time.time() - startTime
print("Elapsed = %.2f (s)" % (elapsedTime))

Elapsed = 4.03 (s)

sampleSizeArray = [int(n) for n in np.logspace(0.0, 5.0)]
expectedConvergence = [1.0 / np.sqrt(n) for n in sampleSizeArray]

title = "Convergence of Monte-Carlo method - problem = %s" % (problem.getName())
graph = ot.Graph(title, "Sample size", "Absolute error", True, "topright")
curve = ot.Cloud(sampleSizeAbsoluteErrorTable, "blue", "fsquare", "")
curve.setLegend("Monte-Carlo")
graph.add(curve)
curve = ot.Curve(sampleSizeArray, expectedConvergence)
curve.setLegend(r"$1/\sqrt{n}$")
graph.add(curve)
graph.setLogScale(ot.GraphImplementation.LOGXY)
graph.setColors(["dodgerblue3", "darkorange1"])
_ = otv.View(graph)

Convergence of Monte-Carlo method - problem = R-S

otv.View.ShowAll()

Total running time of the script: (0 minutes 4.666 seconds)

otbenchmark

Module otbenchmark

Table of Contents

Previous topic

Next topic

This Page

Convergence of Monte-Carlo to estimate the probability in a reliability problem¶