.. DO NOT EDIT.
.. THIS FILE WAS AUTOMATICALLY GENERATED BY SPHINX-GALLERY.
.. TO MAKE CHANGES, EDIT THE SOURCE PYTHON FILE:
.. "auto_data_analysis/sample_analysis/plot_compare_unconditional_conditional_histograms.py"
.. LINE NUMBERS ARE GIVEN BELOW.

.. only:: html

    .. note::
        :class: sphx-glr-download-link-note

        :ref:`Go to the end <sphx_glr_download_auto_data_analysis_sample_analysis_plot_compare_unconditional_conditional_histograms.py>`
        to download the full example code.

.. rst-class:: sphx-glr-example-title

.. _sphx_glr_auto_data_analysis_sample_analysis_plot_compare_unconditional_conditional_histograms.py:


Compare unconditional and conditional histograms
================================================

.. GENERATED FROM PYTHON SOURCE LINES 7-30

In this example, we compare unconditional and conditional histograms for a
simulation.
We consider the :ref:`flooding model<use-case-flood-model>`.
Let :math:`g` be a function which takes four inputs :math:`Q`, :math:`K_s`,
:math:`Z_v` and :math:`Z_m` and returns one output :math:`S`.

We first consider the (unconditional) distribution of the input :math:`Q`.

Let :math:`t` be a given threshold on the output :math:`S`: we consider
the event :math:`S > t`.
Then we consider the conditional distribution of the input :math:`Q` given
that :math:`S > t` that is to say :math:`Q|S > t`.

If these two distributions are significantly different, we conclude that
the input :math:`Q` has an impact on the event :math:`S > t`.

In order to approximate the distribution of the output :math:`S`,
we perform a Monte-Carlo simulation with size 500.
The threshold :math:`t` is chosen as the 90% quantile of the empirical
distribution of :math:`S`.
In this example, the distribution is aproximated by its empirical histogram
(but this could be done with another distribution approximation as well,
such as kernel smoothing for example).

.. GENERATED FROM PYTHON SOURCE LINES 32-38

.. code-block:: Python

    import numpy as np
    from openturns.usecases import flood_model
    import openturns as ot
    import openturns.viewer as viewer









.. GENERATED FROM PYTHON SOURCE LINES 39-40

We use the `FloodModel` data class that contains all the case parameters.

.. GENERATED FROM PYTHON SOURCE LINES 40-43

.. code-block:: Python

    fm = flood_model.FloodModel()









.. GENERATED FROM PYTHON SOURCE LINES 44-46

Create an input sample from the joint `distribution` defined in the data class.
We build an output sample by taking the image by the `model`.

.. GENERATED FROM PYTHON SOURCE LINES 48-52

.. code-block:: Python

    size = 500
    inputSample = fm.distribution.getSample(size)
    inputSample[:5]






.. raw:: html

    <div class="output_subarea output_html rendered_html output_result">
    <table>
      <tr><td></td><th>Q (m3/s)</th><th>Ks</th><th>Zv (m)</th><th>Zm (m)</th><th>B (m)</th><th>L (m)</th><th>Zb (m)</th><th>Hd (m)</th></tr>
      <tr><th>0</th><td>1443.603</td><td>30.15661</td><td>49.11714</td><td>55.59186</td><td>301.3592</td><td>5002.721</td><td>55.58008</td><td>2.081542</td></tr>
      <tr><th>1</th><td>2174.89</td><td>34.6789</td><td>50.76485</td><td>55.87647</td><td>297.3326</td><td>5002.738</td><td>55.51363</td><td>3.596236</td></tr>
      <tr><th>2</th><td>626.1024</td><td>35.75353</td><td>50.0302</td><td>54.66188</td><td>295.8753</td><td>5007.265</td><td>55.6321</td><td>3.345528</td></tr>
      <tr><th>3</th><td>325.8124</td><td>36.66599</td><td>49.02634</td><td>55.36675</td><td>298.1232</td><td>5002.338</td><td>55.56241</td><td>3.514421</td></tr>
      <tr><th>4</th><td>981.3994</td><td>41.10229</td><td>49.39776</td><td>54.84771</td><td>298.9384</td><td>5004.617</td><td>55.5559</td><td>3.610411</td></tr>
    </table>
    </div>
    <br />
    <br />

.. GENERATED FROM PYTHON SOURCE LINES 53-56

.. code-block:: Python

    outputSample = fm.model(inputSample)
    outputSample[:5]






.. raw:: html

    <div class="output_subarea output_html rendered_html output_result">
    <table>
      <tr><td></td><th>H</th><th>S</th><th>C</th></tr>
      <tr><th>0</th><td>2.437721</td><td>-6.106769</td><td>1.010225</td></tr>
      <tr><th>1</th><td>3.102218</td><td>-5.242798</td><td>1.187054</td></tr>
      <tr><th>2</th><td>1.49104</td><td>-7.456382</td><td>0.8211161</td></tr>
      <tr><th>3</th><td>0.89888</td><td>-9.151609</td><td>0.7062944</td></tr>
      <tr><th>4</th><td>1.699544</td><td>-8.069001</td><td>0.7681099</td></tr>
    </table>
    </div>
    <br />
    <br />

.. GENERATED FROM PYTHON SOURCE LINES 57-58

Merge the input and output samples into a single sample.

.. GENERATED FROM PYTHON SOURCE LINES 60-64

.. code-block:: Python

    sample = ot.Sample(inputSample)
    sample.stack(outputSample)
    sample[0:5]






.. raw:: html

    <div class="output_subarea output_html rendered_html output_result">
    <table>
      <tr><td></td><th>Q (m3/s)</th><th>Ks</th><th>Zv (m)</th><th>Zm (m)</th><th>B (m)</th><th>L (m)</th><th>Zb (m)</th><th>Hd (m)</th><th>H</th><th>S</th><th>C</th></tr>
      <tr><th>0</th><td>1443.603</td><td>30.15661</td><td>49.11714</td><td>55.59186</td><td>301.3592</td><td>5002.721</td><td>55.58008</td><td>2.081542</td><td>2.437721</td><td>-6.106769</td><td>1.010225</td></tr>
      <tr><th>1</th><td>2174.89</td><td>34.6789</td><td>50.76485</td><td>55.87647</td><td>297.3326</td><td>5002.738</td><td>55.51363</td><td>3.596236</td><td>3.102218</td><td>-5.242798</td><td>1.187054</td></tr>
      <tr><th>2</th><td>626.1024</td><td>35.75353</td><td>50.0302</td><td>54.66188</td><td>295.8753</td><td>5007.265</td><td>55.6321</td><td>3.345528</td><td>1.49104</td><td>-7.456382</td><td>0.8211161</td></tr>
      <tr><th>3</th><td>325.8124</td><td>36.66599</td><td>49.02634</td><td>55.36675</td><td>298.1232</td><td>5002.338</td><td>55.56241</td><td>3.514421</td><td>0.89888</td><td>-9.151609</td><td>0.7062944</td></tr>
      <tr><th>4</th><td>981.3994</td><td>41.10229</td><td>49.39776</td><td>54.84771</td><td>298.9384</td><td>5004.617</td><td>55.5559</td><td>3.610411</td><td>1.699544</td><td>-8.069001</td><td>0.7681099</td></tr>
    </table>
    </div>
    <br />
    <br />

.. GENERATED FROM PYTHON SOURCE LINES 65-67

Extract the first column of `inputSample` into the sample of the flowrates
:math:`Q`.

.. GENERATED FROM PYTHON SOURCE LINES 69-71

.. code-block:: Python

    sampleQ = inputSample[:, 0]








.. GENERATED FROM PYTHON SOURCE LINES 72-75

The next cell defines a function that computes the conditional sample of
a component given that the a marginal (defined by its index `criteriaComponent`)
exceeds a given threshold, defined by its quantile level.

.. GENERATED FROM PYTHON SOURCE LINES 75-97

.. code-block:: Python



    def computeConditionnedSample(
        sample, alpha=0.9, criteriaComponent=None, selectedComponent=0
    ):
        """
        Return values from the selectedComponent-th component of the sample.
        Selects the values according to the alpha-level quantile of
        the criteriaComponent-th component of the sample.
        """
        dim = sample.getDimension()
        if criteriaComponent is None:
            criteriaComponent = dim - 1
        sortedSample = sample.sortAccordingToAComponent(criteriaComponent)
        quantiles = sortedSample.computeQuantilePerComponent(alpha)
        quantileValue = quantiles[criteriaComponent]
        sortedSampleCriteria = sortedSample[:, criteriaComponent]
        indices = np.where(np.array(sortedSampleCriteria.asPoint()) > quantileValue)[0]
        conditionnedSortedSample = sortedSample[int(indices[0]) :, selectedComponent]
        return conditionnedSortedSample









.. GENERATED FROM PYTHON SOURCE LINES 98-99

Create an histogram for the unconditional flowrates.

.. GENERATED FROM PYTHON SOURCE LINES 101-104

.. code-block:: Python

    numberOfBins = 10
    histogram = ot.HistogramFactory().buildAsHistogram(sampleQ, numberOfBins)








.. GENERATED FROM PYTHON SOURCE LINES 105-106

Extract the sub-sample of the input flowrates `Q` which leads to large values of the output `S`.

.. GENERATED FROM PYTHON SOURCE LINES 108-109

Search the index of the marginal `S` in the columns of the sample.

.. GENERATED FROM PYTHON SOURCE LINES 109-112

.. code-block:: Python

    criteriaComponent = list(sample.getDescription()).index("S")
    criteriaComponent





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

 .. code-block:: none


    9



.. GENERATED FROM PYTHON SOURCE LINES 113-119

.. code-block:: Python

    alpha = 0.9
    selectedComponent = 0
    conditionnedSampleQ = computeConditionnedSample(
        sample, alpha, criteriaComponent, selectedComponent
    )








.. GENERATED FROM PYTHON SOURCE LINES 120-130

We could as well use:

.. code-block::

    # conditionnedHistogram = ot.HistogramFactory().buildAsHistogram(conditionnedSampleQ)

but this creates an histogram with new classes, corresponding
to `conditionnedSampleQ`.
We want to use exactly the same classes as the full sample,
so that the two histograms match.

.. GENERATED FROM PYTHON SOURCE LINES 132-138

.. code-block:: Python

    first = histogram.getFirst()
    width = histogram.getWidth()
    conditionnedHistogram = ot.HistogramFactory().buildAsHistogram(
        conditionnedSampleQ, first, width
    )








.. GENERATED FROM PYTHON SOURCE LINES 139-140

Then creates a graphics with the unconditional and the conditional histograms.

.. GENERATED FROM PYTHON SOURCE LINES 142-151

.. code-block:: Python

    graph = histogram.drawPDF()
    graph.setLegends(["Q"])
    #
    graphConditionnalQ = conditionnedHistogram.drawPDF()
    graphConditionnalQ.setLegends([r"$Q | S > S_{%s}$" % (alpha)])
    graph.add(graphConditionnalQ)
    view = viewer.View(graph)





.. image-sg:: /auto_data_analysis/sample_analysis/images/sphx_glr_plot_compare_unconditional_conditional_histograms_001.svg
   :alt: Q (m3/s) PDF
   :srcset: /auto_data_analysis/sample_analysis/images/sphx_glr_plot_compare_unconditional_conditional_histograms_001.svg
   :class: sphx-glr-single-img





.. GENERATED FROM PYTHON SOURCE LINES 152-153

Show all the graphs.

.. GENERATED FROM PYTHON SOURCE LINES 153-155

.. code-block:: Python

    view.ShowAll()








.. GENERATED FROM PYTHON SOURCE LINES 156-164

We see that the two histograms are very different.
The high values of the input :math:`Q` seem to often lead to a high value of the output :math:`S`.

We could explore this situation further by comparing the unconditional
distribution of :math:`Q` (which is known in this case) with the conditonal
distribution of :math:`Q | S > t`, estimated by kernel smoothing.
This would have the advantage of accuracy, since the kernel smoothing is a
more accurate approximation of a distribution than the histogram.


.. _sphx_glr_download_auto_data_analysis_sample_analysis_plot_compare_unconditional_conditional_histograms.py:

.. only:: html

  .. container:: sphx-glr-footer sphx-glr-footer-example

    .. container:: sphx-glr-download sphx-glr-download-jupyter

      :download:`Download Jupyter notebook: plot_compare_unconditional_conditional_histograms.ipynb <plot_compare_unconditional_conditional_histograms.ipynb>`

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

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

    .. container:: sphx-glr-download sphx-glr-download-zip

      :download:`Download zipped: plot_compare_unconditional_conditional_histograms.zip <plot_compare_unconditional_conditional_histograms.zip>`