.. 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 :ref:`Go to the end ` 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 .. 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-99 .. 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 .. code-block:: none Proba Subset = 0.00035581020000000237 Current coefficient of variation = 0.08805253003879156 .. GENERATED FROM PYTHON SOURCE LINES 100-101 The length of the confidence interval of level :math:`95\%` is: .. GENERATED FROM PYTHON SOURCE LINES 103-106 .. code-block:: default length95 = result.getConfidenceLength() print("Confidence length (0.95) = ", result.getConfidenceLength()) .. rst-class:: sphx-glr-script-out .. code-block:: none Confidence length (0.95) = 0.0001228112975006667 .. GENERATED FROM PYTHON SOURCE LINES 107-108 which enables to build the confidence interval: .. GENERATED FROM PYTHON SOURCE LINES 110-118 .. code-block:: default print( "Confidence interval (0.95) = [", proba - length95 / 2, ", ", proba + length95 / 2, "]", ) .. rst-class:: sphx-glr-script-out .. code-block:: none Confidence interval (0.95) = [ 0.000294404551249669 , 0.0004172158487503357 ] .. GENERATED FROM PYTHON SOURCE LINES 119-120 You can also get the succesive thresholds used by the algorithm: .. GENERATED FROM PYTHON SOURCE LINES 122-125 .. code-block:: default levels = algo.getThresholdPerStep() print("Levels of g = ", levels) .. rst-class:: sphx-glr-script-out .. code-block:: none Levels of g = [56.9472,18.3721,8.502,0] .. GENERATED FROM PYTHON SOURCE LINES 126-132 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 134-138 .. code-block:: default inputSampleSubset = g.getInputHistory() nTotal = inputSampleSubset.getSize() print("Number of evaluations of g = ", nTotal) .. rst-class:: sphx-glr-script-out .. code-block:: none Number of evaluations of g = 40000 .. GENERATED FROM PYTHON SOURCE LINES 139-140 Within each step of the algorithm, a sample of size :math:`N` is created, where: .. GENERATED FROM PYTHON SOURCE LINES 142-145 .. code-block:: default N = algo.getMaximumOuterSampling() * algo.getBlockSize() print("Size of each subset = ", N) .. rst-class:: sphx-glr-script-out .. code-block:: none Size of each subset = 10000 .. GENERATED FROM PYTHON SOURCE LINES 146-147 You can get the number :math:`N_s` of steps with: .. GENERATED FROM PYTHON SOURCE LINES 149-152 .. code-block:: default Ns = algo.getStepsNumber() print("Number of steps= ", Ns) .. rst-class:: sphx-glr-script-out .. code-block:: none Number of steps= 4 .. GENERATED FROM PYTHON SOURCE LINES 153-154 and you can verify that :math:`N_s` is equal to :math:`\frac{nTotal}{N}`: .. GENERATED FROM PYTHON SOURCE LINES 156-158 .. code-block:: default print("nTotal / N = ", int(nTotal / N)) .. rst-class:: sphx-glr-script-out .. code-block:: none nTotal / N = 4 .. GENERATED FROM PYTHON SOURCE LINES 159-160 Now, we can split the initial sample into subset samples of size :math:`N_s`: .. GENERATED FROM PYTHON SOURCE LINES 162-166 .. 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 167-168 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 170-178 .. 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 179-180 Add all the subset samples: .. GENERATED FROM PYTHON SOURCE LINES 182-189 .. 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 190-191 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 193-203 .. 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 204-206 .. 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 .. GENERATED FROM PYTHON SOURCE LINES 207-211 Draw the frontiers only ----------------------- The following graph enables to understand the progresison of the algorithm: .. GENERATED FROM PYTHON SOURCE LINES 213-232 .. 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 .. GENERATED FROM PYTHON SOURCE LINES 233-236 Get all the input and output points that realized the event ----------------------------------------------------------- The following lines are possible only if you have mentioned that you wanted to keep the points that realize the event with the method *algo.setKeepEventSample(True)* .. GENERATED FROM PYTHON SOURCE LINES 238-242 .. code-block:: default inputEventSample = algo.getEventInputSample() outputEventSample = algo.getEventOutputSample() print("Number of event realizations = ", inputEventSample.getSize()) .. rst-class:: sphx-glr-script-out .. code-block:: none Number of event realizations = 3551 .. GENERATED FROM PYTHON SOURCE LINES 243-244 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 246-250 .. 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 251-252 Draw them! They are all in the event space. .. GENERATED FROM PYTHON SOURCE LINES 254-266 .. 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-timing **Total running time of the script:** ( 0 minutes 0.765 seconds) .. _sphx_glr_download_auto_reliability_sensitivity_reliability_plot_subset_sampling.py: .. only:: html .. container:: sphx-glr-footer 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 `