.. 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 <sphx_glr_download_auto_reliability_sensitivity_reliability_plot_subset_sampling.py>`
        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

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:: Python

    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:: Python

    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:: Python

    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:: Python

    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:: Python

    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:: Python

    algo.run()








.. GENERATED FROM PYTHON SOURCE LINES 94-99

.. 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.00040940000000000025
    Current coefficient of variation =  0.08584233339124359




.. 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:: 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.00013776136561433357




.. GENERATED FROM PYTHON SOURCE LINES 107-108

which enables to build the confidence interval:

.. GENERATED FROM PYTHON SOURCE LINES 110-118

.. 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.0003405193171928335 ,  0.000478280682807167 ]




.. GENERATED FROM PYTHON SOURCE LINES 119-120

You can also get the successive thresholds used by the algorithm:

.. GENERATED FROM PYTHON SOURCE LINES 122-125

.. code-block:: Python

    levels = algo.getThresholdPerStep()
    print("Levels of g = ", levels)





.. rst-class:: sphx-glr-script-out

 .. code-block:: none

    Levels of g =  [54.0207,18.0184,7.60975,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:: Python

    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:: 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 =  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:: 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 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:: Python

    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:: Python

    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:: 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("bottomleft")








.. GENERATED FROM PYTHON SOURCE LINES 179-180

Add all the subset samples:

.. GENERATED FROM PYTHON SOURCE LINES 182-189

.. 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 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:: 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 204-206

.. 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 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:: 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("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:: Python

    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 =  4094




.. GENERATED FROM PYTHON SOURCE LINES 243-244

Draw them! They are all in the event space.

.. GENERATED FROM PYTHON SOURCE LINES 246-258

.. code-block:: Python

    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






.. _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 <plot_subset_sampling.ipynb>`

    .. container:: sphx-glr-download sphx-glr-download-python

      :download:`Download Python source code: plot_subset_sampling.py <plot_subset_sampling.py>`