.. DO NOT EDIT. .. THIS FILE WAS AUTOMATICALLY GENERATED BY SPHINX-GALLERY. .. TO MAKE CHANGES, EDIT THE SOURCE PYTHON FILE: .. "auto_reliability_sensitivity/reliability/plot_subset_sampling.py" .. LINE NUMBERS ARE GIVEN BELOW. .. only:: html .. note:: :class: sphx-glr-download-link-note Click :ref:`here ` to download the full example code .. rst-class:: sphx-glr-example-title .. _sphx_glr_auto_reliability_sensitivity_reliability_plot_subset_sampling.py: Subset Sampling =============== .. GENERATED FROM PYTHON SOURCE LINES 6-29 The objective is to evaluate a probability from the Subset sampling technique. We consider the function :math:`g : \mathbb{R}^2 \rightarrow \mathbb{R}` defined by: .. math:: \begin{align*} g(X)= 20-(x_1-x_2)^2-8(x_1+x_2-4)^3 \end{align*} and the input random vector :math:`X = (X_1, X_2)` which follows a Normal distribution with independent components, and identical marginals with 0.25 mean and unit variance: .. math:: \begin{align*} X \sim \mathcal{N}(\mu = [0.25, 0.25], \sigma = [1,1], cov = I_2) \end{align*} We want to evaluate the probability: .. math:: \begin{align*} p = \mathbb{P} \{ g(X) \leq 0 \} \end{align*} .. GENERATED FROM PYTHON SOURCE LINES 32-33 First, import the python modules: .. GENERATED FROM PYTHON SOURCE LINES 35-38 .. code-block:: default import openturns as ot from openturns.viewer import View .. GENERATED FROM PYTHON SOURCE LINES 39-41 Create the probabilistic model :math:`Y = g(X)` ----------------------------------------------- .. GENERATED FROM PYTHON SOURCE LINES 43-44 Create the input random vector :math:`X`: .. GENERATED FROM PYTHON SOURCE LINES 46-48 .. code-block:: default X = ot.RandomVector(ot.Normal([0.25]*2, [1]*2, ot.IdentityMatrix(2))) .. GENERATED FROM PYTHON SOURCE LINES 49-50 Create the function :math:`g`: .. GENERATED FROM PYTHON SOURCE LINES 52-55 .. code-block:: default g = ot.SymbolicFunction(['x1', 'x2'], ['20-(x1-x2)^2-8*(x1+x2-4)^3']) print('function g: ', g) .. rst-class:: sphx-glr-script-out Out: .. code-block:: none function g: [x1,x2]->[20-(x1-x2)^2-8*(x1+x2-4)^3] .. GENERATED FROM PYTHON SOURCE LINES 56-57 In order to be able to get the subset samples used in the algorithm, it is necessary to transform the *SymbolicFunction* into a *MemoizeFunction*: .. GENERATED FROM PYTHON SOURCE LINES 59-61 .. code-block:: default g = ot.MemoizeFunction(g) .. GENERATED FROM PYTHON SOURCE LINES 62-63 Create the output random vector :math:`Y = g(X)`: .. GENERATED FROM PYTHON SOURCE LINES 65-67 .. code-block:: default Y = ot.CompositeRandomVector(g, X) .. GENERATED FROM PYTHON SOURCE LINES 68-70 Create the event :math:`\{ Y = g(X) \leq 0 \}` ---------------------------------------------- .. GENERATED FROM PYTHON SOURCE LINES 72-74 .. code-block:: default myEvent = ot.ThresholdEvent(Y, ot.LessOrEqual(), 0.0) .. GENERATED FROM PYTHON SOURCE LINES 75-77 Evaluate the probability with the subset sampling technique ----------------------------------------------------------- .. GENERATED FROM PYTHON SOURCE LINES 79-81 .. code-block:: default algo = ot.SubsetSampling(myEvent) .. GENERATED FROM PYTHON SOURCE LINES 82-83 In order to get all the inputs and outputs that realize the event, you have to mention it now: .. GENERATED FROM PYTHON SOURCE LINES 85-87 .. code-block:: default algo.setKeepEventSample(True) .. GENERATED FROM PYTHON SOURCE LINES 88-89 Now you can run the algorithm! .. GENERATED FROM PYTHON SOURCE LINES 91-93 .. code-block:: default algo.run() .. GENERATED FROM PYTHON SOURCE LINES 94-100 .. code-block:: default result = algo.getResult() proba = result.getProbabilityEstimate() print('Proba Subset = ', proba) print('Current coefficient of variation = ', result.getCoefficientOfVariation()) .. rst-class:: sphx-glr-script-out Out: .. code-block:: none Proba Subset = 0.00042720000000000176 Current coefficient of variation = 0.08755460995974475 .. GENERATED FROM PYTHON SOURCE LINES 101-102 The length of the confidence interval of level :math:`95\%` is: .. GENERATED FROM PYTHON SOURCE LINES 104-107 .. code-block:: default length95 = result.getConfidenceLength() print('Confidence length (0.95) = ', result.getConfidenceLength()) .. rst-class:: sphx-glr-script-out Out: .. code-block:: none Confidence length (0.95) = 0.00014661835695300627 .. GENERATED FROM PYTHON SOURCE LINES 108-109 which enables to build the confidence interval: .. GENERATED FROM PYTHON SOURCE LINES 111-114 .. code-block:: default print('Confidence intervalle (0.95) = [', proba - length95/2, ', ', proba + length95/2, ']') .. rst-class:: sphx-glr-script-out Out: .. code-block:: none Confidence intervalle (0.95) = [ 0.00035389082152349863 , 0.000500509178476505 ] .. GENERATED FROM PYTHON SOURCE LINES 115-116 You can also get the succesive thresholds used by the algorithm: .. GENERATED FROM PYTHON SOURCE LINES 118-121 .. code-block:: default levels = algo.getThresholdPerStep() print('Levels of g = ', levels) .. rst-class:: sphx-glr-script-out Out: .. code-block:: none Levels of g = [57.5557,18.1503,7.37784,0] .. GENERATED FROM PYTHON SOURCE LINES 122-128 Draw the subset samples used by the algorithm --------------------------------------------- The following manipulations are possible onfly if you have created a *MemoizeFunction* that enables to store all the inputs and output of the function :math:`g`. Get all the inputs where :math:`g` were evaluated: .. GENERATED FROM PYTHON SOURCE LINES 130-134 .. code-block:: default inputSampleSubset = g.getInputHistory() nTotal = inputSampleSubset.getSize() print('Number of evaluations of g = ', nTotal) .. rst-class:: sphx-glr-script-out Out: .. code-block:: none Number of evaluations of g = 40000 .. GENERATED FROM PYTHON SOURCE LINES 135-136 Within each step of the algorithm, a sample of size :math:`N` is created, where: .. GENERATED FROM PYTHON SOURCE LINES 138-141 .. code-block:: default N = algo.getMaximumOuterSampling()*algo.getBlockSize() print('Size of each subset = ', N) .. rst-class:: sphx-glr-script-out Out: .. code-block:: none Size of each subset = 10000 .. GENERATED FROM PYTHON SOURCE LINES 142-143 You can get the number :math:`N_s` of steps with: .. GENERATED FROM PYTHON SOURCE LINES 145-148 .. code-block:: default Ns = algo.getStepsNumber() print('Number of steps= ', Ns) .. rst-class:: sphx-glr-script-out Out: .. code-block:: none Number of steps= 4 .. GENERATED FROM PYTHON SOURCE LINES 149-150 and you can verify that :math:`N_s` is equal to :math:`\frac{nTotal}{N}`: .. GENERATED FROM PYTHON SOURCE LINES 152-154 .. code-block:: default print('nTotal / N = ', int(nTotal / N)) .. rst-class:: sphx-glr-script-out Out: .. code-block:: none nTotal / N = 4 .. GENERATED FROM PYTHON SOURCE LINES 155-156 Now, we can split the initial sample into subset samples of size :math:`N_s`: .. GENERATED FROM PYTHON SOURCE LINES 158-162 .. code-block:: default list_subSamples = list() for i in range(Ns): list_subSamples.append(inputSampleSubset[i*N:i*N + N]) .. GENERATED FROM PYTHON SOURCE LINES 163-164 The following graph draws each subset sample and the frontier :math:`g(x_1, x_2) = l_i` where :math:`l_i` is the threshold at the step :math:`i`: .. GENERATED FROM PYTHON SOURCE LINES 166-174 .. code-block:: default graph = ot.Graph() graph.setAxes(True) graph.setGrid(True) graph.setTitle('Subset sampling: samples') graph.setXTitle(r'$x_1$') graph.setYTitle(r'$x_2$') graph.setLegendPosition('bottomleft') .. GENERATED FROM PYTHON SOURCE LINES 175-176 Add all the subset samples: .. GENERATED FROM PYTHON SOURCE LINES 178-185 .. code-block:: default for i in range(Ns): cloud = ot.Cloud(list_subSamples[i]) # cloud.setPointStyle("dot") graph.add(cloud) col = ot.Drawable().BuildDefaultPalette(Ns) graph.setColors(col) .. GENERATED FROM PYTHON SOURCE LINES 186-187 Add the frontiers :math:`g(x_1, x_2) = l_i` where :math:`l_i` is the threshold at the step :math:`i`: .. GENERATED FROM PYTHON SOURCE LINES 189-199 .. code-block:: default gIsoLines = g.draw([-3]*2, [5]*2, [128]*2) dr = gIsoLines.getDrawable(0) for i in range(levels.getSize()): dr.setLevels([levels[i]]) dr.setLineStyle('solid') dr.setLegend(r'$g(X) = $' + str(round(levels[i], 2))) dr.setLineWidth(3) dr.setColor(col[i]) graph.add(dr) .. GENERATED FROM PYTHON SOURCE LINES 200-202 .. code-block:: default View(graph) .. image-sg:: /auto_reliability_sensitivity/reliability/images/sphx_glr_plot_subset_sampling_001.png :alt: Subset sampling: samples :srcset: /auto_reliability_sensitivity/reliability/images/sphx_glr_plot_subset_sampling_001.png :class: sphx-glr-single-img .. rst-class:: sphx-glr-script-out Out: .. code-block:: none .. GENERATED FROM PYTHON SOURCE LINES 203-207 Draw the frontiers only ----------------------- The following graph enables to understand the progresison of the algorithm: .. GENERATED FROM PYTHON SOURCE LINES 209-228 .. code-block:: default graph = ot.Graph() graph.setAxes(True) graph.setGrid(True) dr = gIsoLines.getDrawable(0) for i in range(levels.getSize()): dr.setLevels([levels[i]]) dr.setLineStyle('solid') dr.setLegend(r'$g(X) = $' + str(round(levels[i], 2))) dr.setLineWidth(3) graph.add(dr) graph.setColors(col) graph.setLegendPosition('bottomleft') graph.setTitle('Subset sampling: thresholds') graph.setXTitle(r'$x_1$') graph.setYTitle(r'$x_2$') View(graph) .. image-sg:: /auto_reliability_sensitivity/reliability/images/sphx_glr_plot_subset_sampling_002.png :alt: Subset sampling: thresholds :srcset: /auto_reliability_sensitivity/reliability/images/sphx_glr_plot_subset_sampling_002.png :class: sphx-glr-single-img .. rst-class:: sphx-glr-script-out Out: .. code-block:: none .. GENERATED FROM PYTHON SOURCE LINES 229-232 Get all the input and output points that realized the event ----------------------------------------------------------- The following lines are possible only if you have mentionned that you wanted to keep the points that realize the event with the method *algo.setKeepEventSample(True)* .. GENERATED FROM PYTHON SOURCE LINES 234-238 .. code-block:: default inputEventSample = algo.getEventInputSample() outputEventSample = algo.getEventOutputSample() print('Number of event realizations = ', inputEventSample.getSize()) .. rst-class:: sphx-glr-script-out Out: .. code-block:: none Number of event realizations = 4272 .. GENERATED FROM PYTHON SOURCE LINES 239-240 Here we have to avoid a bug of the version 1.15 because *getEventInputSample()* gives the sample in the stadrad space: we have to push it backward to the physical space. .. GENERATED FROM PYTHON SOURCE LINES 242-246 .. code-block:: default dist = ot.Normal([0.25]*2, [1]*2, ot.IdentityMatrix(2)) transformFunc = dist.getInverseIsoProbabilisticTransformation() inputEventSample = transformFunc(inputEventSample) .. GENERATED FROM PYTHON SOURCE LINES 247-248 Draw them! They are all in the event space. .. GENERATED FROM PYTHON SOURCE LINES 250-262 .. code-block:: default graph = ot.Graph() graph.setAxes(True) graph.setGrid(True) cloud = ot.Cloud(inputEventSample) cloud.setPointStyle('dot') graph.add(cloud) gIsoLines = g.draw([-3]*2, [5]*2, [1000]*2) dr = gIsoLines.getDrawable(0) dr.setLevels([0.0]) dr.setColor('red') graph.add(dr) View(graph) .. image-sg:: /auto_reliability_sensitivity/reliability/images/sphx_glr_plot_subset_sampling_003.png :alt: plot subset sampling :srcset: /auto_reliability_sensitivity/reliability/images/sphx_glr_plot_subset_sampling_003.png :class: sphx-glr-single-img .. rst-class:: sphx-glr-script-out Out: .. code-block:: none .. rst-class:: sphx-glr-timing **Total running time of the script:** ( 0 minutes 1.112 seconds) .. _sphx_glr_download_auto_reliability_sensitivity_reliability_plot_subset_sampling.py: .. only :: html .. container:: sphx-glr-footer :class: sphx-glr-footer-example .. container:: sphx-glr-download sphx-glr-download-python :download:`Download Python source code: plot_subset_sampling.py ` .. container:: sphx-glr-download sphx-glr-download-jupyter :download:`Download Jupyter notebook: plot_subset_sampling.ipynb ` .. only:: html .. rst-class:: sphx-glr-signature `Gallery generated by Sphinx-Gallery `_