.. DO NOT EDIT. .. THIS FILE WAS AUTOMATICALLY GENERATED BY SPHINX-GALLERY. .. TO MAKE CHANGES, EDIT THE SOURCE PYTHON FILE: .. "auto_reliability_sensitivity/reliability/plot_nais.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_nais.py: Non parametric Adaptive Importance Sampling (NAIS) ================================================== .. GENERATED FROM PYTHON SOURCE LINES 6-33 The objective is to evaluate a probability from the Non parametric Adaptive Importance Sampling (NAIS) technique. We consider the four-branch function :math:`g : \mathbb{R}^2 \rightarrow \mathbb{R}` defined by: .. math:: \begin{align*} g(\vect{X}) = \min \begin{pmatrix}5+0.1(x_1-x_2)^2-\frac{(x_1+x_2)}{\sqrt{2}}\\ 5+0.1(x_1-x_2)^2+\frac{(x_1+x_2)}{\sqrt{2}}\\ (x_1-x_2)+ \frac{9}{\sqrt{2}}\\ (x_2-x_1)+ \frac{9}{\sqrt{2}} \end{pmatrix} \end{align*} and the input random vector :math:`\vect{X} = (X_1, X_2)` which follows the standard 2-dimensional Normal distribution: .. math:: \begin{align*} \vect{X} \sim \mathcal{N}(\mu = [0, 0], \sigma = [1,1], corr = \mat{I}_2) \end{align*} We want to evaluate the probability: .. math:: \begin{align*} p = \mathbb{P} ( g(\vect{X}) \leq 0 ) \end{align*} .. GENERATED FROM PYTHON SOURCE LINES 36-37 First, import the python modules: .. GENERATED FROM PYTHON SOURCE LINES 39-43 .. code-block:: Python import openturns as ot from openturns.viewer import View import math .. GENERATED FROM PYTHON SOURCE LINES 44-46 Create the probabilistic model :math:`Y = g(\vect{X})` ------------------------------------------------------ .. GENERATED FROM PYTHON SOURCE LINES 48-49 Create the input random vector :math:`\vect{X}`: .. GENERATED FROM PYTHON SOURCE LINES 51-53 .. code-block:: Python X = ot.RandomVector(ot.Normal(2)) .. GENERATED FROM PYTHON SOURCE LINES 54-55 Create the function :math:`g` from a :class:`~openturns.PythonFunction`: .. GENERATED FROM PYTHON SOURCE LINES 57-73 .. code-block:: Python def fourBranch(x): x1 = x[0] x2 = x[1] g1 = 5 + 0.1 * (x1 - x2) ** 2 - (x1 + x2) / math.sqrt(2) g2 = 5 + 0.1 * (x1 - x2) ** 2 + (x1 + x2) / math.sqrt(2) g3 = (x1 - x2) + 9 / math.sqrt(2) g4 = (x2 - x1) + 9 / math.sqrt(2) return [min((g1, g2, g3, g4))] g = ot.PythonFunction(2, 1, fourBranch) .. GENERATED FROM PYTHON SOURCE LINES 74-75 Draw the function :math:`g` to help to understand the shape of the limit state function: .. GENERATED FROM PYTHON SOURCE LINES 77-83 .. code-block:: Python graph = ot.Graph("Four Branch function", "x1", "x2", True, "upper right") drawfunction = g.draw([-8] * 2, [8] * 2, [100] * 2) graph.add(drawfunction) view = View(graph) .. image-sg:: /auto_reliability_sensitivity/reliability/images/sphx_glr_plot_nais_001.png :alt: Four Branch function :srcset: /auto_reliability_sensitivity/reliability/images/sphx_glr_plot_nais_001.png :class: sphx-glr-single-img .. GENERATED FROM PYTHON SOURCE LINES 84-85 In order to be able to get the NAIS samples used in the algorithm, it is necessary to transform the :class:`~openturns.PythonFunction` into a :class:`~openturns.MemoizeFunction`: .. GENERATED FROM PYTHON SOURCE LINES 87-89 .. code-block:: Python g = ot.MemoizeFunction(g) .. GENERATED FROM PYTHON SOURCE LINES 90-91 Create the output random vector :math:`Y = g(\vect{X})`: .. GENERATED FROM PYTHON SOURCE LINES 93-95 .. code-block:: Python Y = ot.CompositeRandomVector(g, X) .. GENERATED FROM PYTHON SOURCE LINES 96-98 Create the event :math:`\{ Y = g(\vect{X}) \leq 0 \}` ----------------------------------------------------- .. GENERATED FROM PYTHON SOURCE LINES 100-103 .. code-block:: Python threshold = 0.0 myEvent = ot.ThresholdEvent(Y, ot.Less(), threshold) .. GENERATED FROM PYTHON SOURCE LINES 104-106 Evaluate the probability with the NAIS technique ------------------------------------------------ .. GENERATED FROM PYTHON SOURCE LINES 108-111 .. code-block:: Python quantileLevel = 0.1 algo = ot.NAIS(myEvent, quantileLevel) .. GENERATED FROM PYTHON SOURCE LINES 112-113 Now you can run the algorithm. .. GENERATED FROM PYTHON SOURCE LINES 115-121 .. code-block:: Python algo.run() result = algo.getResult() proba = result.getProbabilityEstimate() print("Proba NAIS = ", proba) print("Current coefficient of variation = ", result.getCoefficientOfVariation()) .. rst-class:: sphx-glr-script-out .. code-block:: none Proba NAIS = 8.001913221123705e-06 Current coefficient of variation = 0.09906223602270714 .. GENERATED FROM PYTHON SOURCE LINES 122-123 The length of the confidence interval of level :math:`95\%` is: .. GENERATED FROM PYTHON SOURCE LINES 125-128 .. 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) = 3.1072775732814027e-06 .. GENERATED FROM PYTHON SOURCE LINES 129-130 which enables to build the confidence interval: .. GENERATED FROM PYTHON SOURCE LINES 132-140 .. 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) = [ 6.4482744344830036e-06 , 9.555552007764406e-06 ] .. GENERATED FROM PYTHON SOURCE LINES 141-147 Draw the NAIS samples used by the algorithm ------------------------------------------- The following manipulations are possible only if you have created a :class:`~openturns.MemoizeFunction` that enables to store all the inputs and outputs of the function :math:`g`. Get all the inputs and outputs where :math:`g` were evaluated: .. GENERATED FROM PYTHON SOURCE LINES 149-154 .. code-block:: Python inputNAIS = g.getInputHistory() outputNAIS = g.getOutputHistory() nTotal = inputNAIS.getSize() print("Number of evaluations of g = ", nTotal) .. rst-class:: sphx-glr-script-out .. code-block:: none Number of evaluations of g = 4000 .. GENERATED FROM PYTHON SOURCE LINES 155-156 Within each step of the algorithm, a sample of size :math:`N` is created, where: .. GENERATED FROM PYTHON SOURCE LINES 158-161 .. code-block:: Python N = algo.getMaximumOuterSampling() * algo.getBlockSize() print("Size of each subset = ", N) .. rst-class:: sphx-glr-script-out .. code-block:: none Size of each subset = 1000 .. GENERATED FROM PYTHON SOURCE LINES 162-163 You can get the number :math:`N_s` of steps with: .. GENERATED FROM PYTHON SOURCE LINES 165-168 .. code-block:: Python Ns = int(nTotal / N) print("Number of steps = ", Ns) .. rst-class:: sphx-glr-script-out .. code-block:: none Number of steps = 4 .. GENERATED FROM PYTHON SOURCE LINES 169-170 Now, we can split the initial sample into NAIS samples of size :math:`N_s`: .. GENERATED FROM PYTHON SOURCE LINES 172-178 .. code-block:: Python listNAISSamples = list() listOutputNAISSamples = list() for i in range(Ns): listNAISSamples.append(inputNAIS[i * N: i * N + N]) listOutputNAISSamples.append(outputNAIS[i * N: i * N + N]) .. GENERATED FROM PYTHON SOURCE LINES 179-180 And get all the levels defining the intermediate and final thresholds given by the empirical quantiles of each NAIS output sample: .. GENERATED FROM PYTHON SOURCE LINES 182-187 .. code-block:: Python levels = [] for i in range(Ns - 1): levels.append(listOutputNAISSamples[i].computeQuantile(quantileLevel)[0]) levels.append(threshold) .. GENERATED FROM PYTHON SOURCE LINES 188-189 The following graph draws each NAIS 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 191-194 .. code-block:: Python graph = ot.Graph("NAIS samples", "x1", "x2", True, "lower left") graph.setGrid(True) .. GENERATED FROM PYTHON SOURCE LINES 195-196 Add all the NAIS samples: .. GENERATED FROM PYTHON SOURCE LINES 198-204 .. code-block:: Python for i in range(Ns): cloud = ot.Cloud(listNAISSamples[i]) graph.add(cloud) col = ot.Drawable().BuildDefaultPalette(Ns) graph.setColors(col) .. GENERATED FROM PYTHON SOURCE LINES 205-206 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 208-221 .. code-block:: Python gIsoLines = g.draw([-8] * 2, [8] * 2, [128] * 2) dr = gIsoLines.getDrawable(2) for i, lv in enumerate(levels): dr.setLevels([lv]) dr.setLineStyle("solid") dr.setLegend(r"$g(X) = $" + str(round(lv, 2))) dr.setLineWidth(3) dr.setColor(col[i]) graph.add(dr) # sphinx_gallery_thumbnail_number = 2 view = View(graph) .. image-sg:: /auto_reliability_sensitivity/reliability/images/sphx_glr_plot_nais_002.png :alt: NAIS samples :srcset: /auto_reliability_sensitivity/reliability/images/sphx_glr_plot_nais_002.png :class: sphx-glr-single-img .. GENERATED FROM PYTHON SOURCE LINES 222-226 Draw the frontiers only ----------------------- The following graph enables to understand the progression of the algorithm from the mean value of the initial distribution to the limit state function: .. GENERATED FROM PYTHON SOURCE LINES 228-240 .. code-block:: Python graph = ot.Graph("NAIS thresholds", "x1", "x2", True, "lower left") graph.setGrid(True) dr = gIsoLines.getDrawable(0) for i, lv in enumerate(levels): dr.setLevels([lv]) dr.setLineStyle("solid") dr.setLegend(r"$g(X) = $" + str(round(lv, 2))) dr.setLineWidth(3) graph.add(dr) graph.setColors(col) view = View(graph) .. image-sg:: /auto_reliability_sensitivity/reliability/images/sphx_glr_plot_nais_003.png :alt: NAIS thresholds :srcset: /auto_reliability_sensitivity/reliability/images/sphx_glr_plot_nais_003.png :class: sphx-glr-single-img .. _sphx_glr_download_auto_reliability_sensitivity_reliability_plot_nais.py: .. only:: html .. container:: sphx-glr-footer sphx-glr-footer-example .. container:: sphx-glr-download sphx-glr-download-jupyter :download:`Download Jupyter notebook: plot_nais.ipynb ` .. container:: sphx-glr-download sphx-glr-download-python :download:`Download Python source code: plot_nais.py `