.. 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:: Python 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:: Python 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:: Python 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 Create the output random vector :math:`Y = g(X)`: .. GENERATED FROM PYTHON SOURCE LINES 59-61 .. code-block:: Python Y = ot.CompositeRandomVector(g, X) .. GENERATED FROM PYTHON SOURCE LINES 62-64 Create the event :math:`\{ Y = g(X) \leq 0 \}` ---------------------------------------------- .. GENERATED FROM PYTHON SOURCE LINES 66-68 .. code-block:: Python myEvent = ot.ThresholdEvent(Y, ot.LessOrEqual(), 0.0) .. GENERATED FROM PYTHON SOURCE LINES 69-71 Evaluate the probability with the subset sampling technique ----------------------------------------------------------- .. GENERATED FROM PYTHON SOURCE LINES 73-75 .. code-block:: Python algo = ot.SubsetSampling(myEvent) .. GENERATED FROM PYTHON SOURCE LINES 76-77 In order to get all the inputs and outputs that realize the event, you have to mention it now: .. GENERATED FROM PYTHON SOURCE LINES 79-81 .. code-block:: Python algo.setKeepSample(True) .. GENERATED FROM PYTHON SOURCE LINES 82-83 Now you can run the algorithm! .. GENERATED FROM PYTHON SOURCE LINES 85-87 .. code-block:: Python algo.run() .. GENERATED FROM PYTHON SOURCE LINES 88-93 .. code-block:: Python 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.0004055999999999982 Current coefficient of variation = 0.09031574730348539 .. GENERATED FROM PYTHON SOURCE LINES 94-95 The length of the confidence interval of level :math:`95\%` is: .. GENERATED FROM PYTHON SOURCE LINES 97-100 .. code-block:: Python length95 = result.getConfidenceLength() print("Confidence length (0.95) = ", result.getConfidenceLength()) .. rst-class:: sphx-glr-script-out .. code-block:: none Confidence length (0.95) = 0.00014359506441517935 .. GENERATED FROM PYTHON SOURCE LINES 101-102 which enables to build the confidence interval: .. GENERATED FROM PYTHON SOURCE LINES 104-112 .. code-block:: Python 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.00033380246779240856 , 0.00047739753220758785 ] .. GENERATED FROM PYTHON SOURCE LINES 113-114 You can also get the successive thresholds used by the algorithm: .. GENERATED FROM PYTHON SOURCE LINES 116-119 .. code-block:: Python levels = algo.getThresholdPerStep() print("Levels of g = ", levels) .. rst-class:: sphx-glr-script-out .. code-block:: none Levels of g = [55.0962,18.3439,7.97654,0] .. GENERATED FROM PYTHON SOURCE LINES 120-122 Draw the subset samples used by the algorithm --------------------------------------------- .. GENERATED FROM PYTHON SOURCE LINES 124-125 You can get the number :math:`N_s` of steps with: .. GENERATED FROM PYTHON SOURCE LINES 125-128 .. code-block:: Python 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 129-130 Get all the inputs where :math:`g` was evaluated at each step .. GENERATED FROM PYTHON SOURCE LINES 130-134 .. code-block:: Python list_subSamples = list() for step in range(Ns): list_subSamples.append(algo.getInputSample(step)) .. GENERATED FROM PYTHON SOURCE LINES 135-136 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 138-146 .. code-block:: Python 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("lower left") .. GENERATED FROM PYTHON SOURCE LINES 147-148 Add all the subset samples: .. GENERATED FROM PYTHON SOURCE LINES 150-157 .. code-block:: Python 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 158-159 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 161-171 .. code-block:: Python 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 172-174 .. code-block:: Python _ = 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 175-179 Draw the frontiers only ----------------------- The following graph enables to understand the progresison of the algorithm: .. GENERATED FROM PYTHON SOURCE LINES 181-200 .. code-block:: Python 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("lower left") 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 201-204 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 samples with the method *algo.setKeepSample(True)* .. GENERATED FROM PYTHON SOURCE LINES 206-212 .. code-block:: Python select = ot.SubsetSampling.EVENT1 # points that realize the event step = Ns - 1 # get the working sample from last iteration inputEventSample = algo.getInputSample(step, select) outputEventSample = algo.getOutputSample(step, select) print("Number of event realizations = ", inputEventSample.getSize()) .. rst-class:: sphx-glr-script-out .. code-block:: none Number of event realizations = 4056 .. GENERATED FROM PYTHON SOURCE LINES 213-214 Draw them! They are all in the event space. .. GENERATED FROM PYTHON SOURCE LINES 216-230 .. code-block:: Python graph = ot.Graph() graph.setAxes(True) graph.setGrid(True) cloud = ot.Cloud(inputEventSample) cloud.setPointStyle("bullet") 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) View.ShowAll() .. 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 .. _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-jupyter :download:`Download Jupyter notebook: plot_subset_sampling.ipynb ` .. container:: sphx-glr-download sphx-glr-download-python :download:`Download Python source code: plot_subset_sampling.py `