.. DO NOT EDIT.
.. THIS FILE WAS AUTOMATICALLY GENERATED BY SPHINX-GALLERY.
.. TO MAKE CHANGES, EDIT THE SOURCE PYTHON FILE:
.. "auto_reliability_sensitivity/sensitivity_analysis/plot_hsic_estimators_ishigami.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_sensitivity_analysis_plot_hsic_estimators_ishigami.py>`
        to download the full example code

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

.. _sphx_glr_auto_reliability_sensitivity_sensitivity_analysis_plot_hsic_estimators_ishigami.py:


The HSIC sensitivity indices: the Ishigami model
================================================

.. GENERATED FROM PYTHON SOURCE LINES 5-10

.. code-block:: Python

    import openturns as ot
    import openturns.viewer as otv
    from openturns.usecases import ishigami_function









.. GENERATED FROM PYTHON SOURCE LINES 11-14

This example is a brief overview of the HSIC sensitivity indices classes and how to call them.
HSIC estimators rely on a reproducing kernel of a Hilbert space. We can use them to compute sensitivity
indices. We present the methods on the :ref:`Ishigami function<use-case-ishigami>`.

.. GENERATED FROM PYTHON SOURCE LINES 18-22

Definition of the model
-----------------------

We load the model from the usecases module.

.. GENERATED FROM PYTHON SOURCE LINES 22-24

.. code-block:: Python

    im = ishigami_function.IshigamiModel()








.. GENERATED FROM PYTHON SOURCE LINES 25-26

We generate an input sample of size 100 (and dimension 3).

.. GENERATED FROM PYTHON SOURCE LINES 26-29

.. code-block:: Python

    size = 100
    X = im.distributionX.getSample(size)








.. GENERATED FROM PYTHON SOURCE LINES 30-31

We compute the output by applying the Ishigami model to the input sample.

.. GENERATED FROM PYTHON SOURCE LINES 31-34

.. code-block:: Python

    Y = im.model(X)









.. GENERATED FROM PYTHON SOURCE LINES 35-45

Setting the covariance models
-----------------------------

The HSIC algorithms use reproducing kernels defined on Hilbert spaces to estimate independence.
For each input variable we choose a covariance kernel.
Here we choose a :class:`~openturns.SquaredExponential`
kernel for all input variables.

They are all stored in a list of :math:`d+1` covariance kernels where :math:`d` is the number of
input variables. The remaining one is for the output variable.

.. GENERATED FROM PYTHON SOURCE LINES 45-47

.. code-block:: Python

    covarianceModelCollection = []








.. GENERATED FROM PYTHON SOURCE LINES 48-54

.. code-block:: Python

    for i in range(3):
        Xi = X.getMarginal(i)
        inputCovariance = ot.SquaredExponential(1)
        inputCovariance.setScale(Xi.computeStandardDeviation())
        covarianceModelCollection.append(inputCovariance)








.. GENERATED FROM PYTHON SOURCE LINES 55-56

Likewise we define a covariance kernel associated to the output variable.

.. GENERATED FROM PYTHON SOURCE LINES 56-60

.. code-block:: Python

    outputCovariance = ot.SquaredExponential(1)
    outputCovariance.setScale(Y.computeStandardDeviation())
    covarianceModelCollection.append(outputCovariance)








.. GENERATED FROM PYTHON SOURCE LINES 61-66

The Global HSIC estimator
-------------------------

In this paragraph, we perform the analysis on the raw data: that is
the global HSIC estimator.

.. GENERATED FROM PYTHON SOURCE LINES 68-80

Choosing an estimator
^^^^^^^^^^^^^^^^^^^^^

After having defined the covariance kernels one has to select an
appropriate estimator for the computations.

Two estimators are proposed:

- an unbiased estimator through the :class:`~openturns.HSICUStat` class
- a biased, but asymptotically unbiased, estimator through the :class:`~openturns.HSICVStat` class

Beware that the conditional analysis used later cannot be performed with the unbiased estimator.

.. GENERATED FROM PYTHON SOURCE LINES 80-83

.. code-block:: Python

    estimatorType = ot.HSICUStat()









.. GENERATED FROM PYTHON SOURCE LINES 84-85

We now build the HSIC estimator:

.. GENERATED FROM PYTHON SOURCE LINES 85-89

.. code-block:: Python

    globHSIC = ot.HSICEstimatorGlobalSensitivity(
        covarianceModelCollection, X, Y, estimatorType
    )








.. GENERATED FROM PYTHON SOURCE LINES 90-91

We get the R2-HSIC indices:

.. GENERATED FROM PYTHON SOURCE LINES 91-95

.. code-block:: Python

    R2HSICIndices = globHSIC.getR2HSICIndices()
    print("\n Global HSIC analysis")
    print("R2-HSIC Indices: ", R2HSICIndices)





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

 .. code-block:: none


     Global HSIC analysis
    R2-HSIC Indices:  [0.249305,-0.00429972,0.0437032]




.. GENERATED FROM PYTHON SOURCE LINES 96-97

and the HSIC indices:

.. GENERATED FROM PYTHON SOURCE LINES 97-100

.. code-block:: Python

    HSICIndices = globHSIC.getHSICIndices()
    print("HSIC Indices: ", HSICIndices)





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

 .. code-block:: none

    HSIC Indices:  [0.0204961,-0.000366135,0.00366669]




.. GENERATED FROM PYTHON SOURCE LINES 101-102

The p-value by permutation.

.. GENERATED FROM PYTHON SOURCE LINES 102-105

.. code-block:: Python

    pvperm = globHSIC.getPValuesPermutation()
    print("p-value (permutation): ", pvperm)





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

 .. code-block:: none

    p-value (permutation):  [0,0.50495,0.00990099]




.. GENERATED FROM PYTHON SOURCE LINES 106-107

We have an asymptotic estimate of the value for this estimator.

.. GENERATED FROM PYTHON SOURCE LINES 107-110

.. code-block:: Python

    pvas = globHSIC.getPValuesAsymptotic()
    print("p-value (asymptotic): ", pvas)





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

 .. code-block:: none

    p-value (asymptotic):  [4.62161e-12,0.553716,0.0159241]




.. GENERATED FROM PYTHON SOURCE LINES 111-112

We vizualise the results.

.. GENERATED FROM PYTHON SOURCE LINES 112-125

.. code-block:: Python

    graph1 = globHSIC.drawHSICIndices()
    view1 = otv.View(graph1)

    graph2 = globHSIC.drawPValuesAsymptotic()
    view2 = otv.View(graph2)

    graph3 = globHSIC.drawR2HSICIndices()
    view3 = otv.View(graph3)

    graph4 = globHSIC.drawPValuesPermutation()
    view4 = otv.View(graph4)





.. rst-class:: sphx-glr-horizontal


    *

      .. image-sg:: /auto_reliability_sensitivity/sensitivity_analysis/images/sphx_glr_plot_hsic_estimators_ishigami_001.png
         :alt: HSIC indices
         :srcset: /auto_reliability_sensitivity/sensitivity_analysis/images/sphx_glr_plot_hsic_estimators_ishigami_001.png
         :class: sphx-glr-multi-img

    *

      .. image-sg:: /auto_reliability_sensitivity/sensitivity_analysis/images/sphx_glr_plot_hsic_estimators_ishigami_002.png
         :alt: Asymptotic p-values
         :srcset: /auto_reliability_sensitivity/sensitivity_analysis/images/sphx_glr_plot_hsic_estimators_ishigami_002.png
         :class: sphx-glr-multi-img

    *

      .. image-sg:: /auto_reliability_sensitivity/sensitivity_analysis/images/sphx_glr_plot_hsic_estimators_ishigami_003.png
         :alt: R2-HSIC indices
         :srcset: /auto_reliability_sensitivity/sensitivity_analysis/images/sphx_glr_plot_hsic_estimators_ishigami_003.png
         :class: sphx-glr-multi-img

    *

      .. image-sg:: /auto_reliability_sensitivity/sensitivity_analysis/images/sphx_glr_plot_hsic_estimators_ishigami_004.png
         :alt: p-values by permutation
         :srcset: /auto_reliability_sensitivity/sensitivity_analysis/images/sphx_glr_plot_hsic_estimators_ishigami_004.png
         :class: sphx-glr-multi-img





.. GENERATED FROM PYTHON SOURCE LINES 126-131

The Target HSIC estimator
-------------------------

We now perform the target analysis which consists in using a filter function over the
output.

.. GENERATED FROM PYTHON SOURCE LINES 134-140

Defining a filter function
^^^^^^^^^^^^^^^^^^^^^^^^^^

We define a filter function on the output variable for the target
analysis. In this example we use the function :math:`\exp{(-d/s)}` where :math:`d` is the distance
to a well-chosen interval.

.. GENERATED FROM PYTHON SOURCE LINES 142-143

We first define a critical domain: in our case that is the :math:`[5,+\infty[` interval.

.. GENERATED FROM PYTHON SOURCE LINES 143-145

.. code-block:: Python

    criticalDomain = ot.Interval(5, float("inf"))








.. GENERATED FROM PYTHON SOURCE LINES 146-148

We have access to the distance to this  domain thanks to the
:class:`~openturns.DistanceToDomainFunction` class.

.. GENERATED FROM PYTHON SOURCE LINES 148-150

.. code-block:: Python

    dist2criticalDomain = ot.DistanceToDomainFunction(criticalDomain)








.. GENERATED FROM PYTHON SOURCE LINES 151-152

We define the parameters in our function from the output sample

.. GENERATED FROM PYTHON SOURCE LINES 152-154

.. code-block:: Python

    s = 0.1 * Y.computeStandardDeviation()[0]








.. GENERATED FROM PYTHON SOURCE LINES 155-157

We now define our filter function by composition of the parametrized function and
the distance function.

.. GENERATED FROM PYTHON SOURCE LINES 157-162

.. code-block:: Python

    f = ot.SymbolicFunction(["x", "s"], ["exp(-x/s)"])
    phi = ot.ParametricFunction(f, [1], [s])
    filterFunction = ot.ComposedFunction(phi, dist2criticalDomain)









.. GENERATED FROM PYTHON SOURCE LINES 163-164

We choose an unbiased estimator

.. GENERATED FROM PYTHON SOURCE LINES 164-166

.. code-block:: Python

    estimatorType = ot.HSICUStat()








.. GENERATED FROM PYTHON SOURCE LINES 167-168

and build the HSIC estimator

.. GENERATED FROM PYTHON SOURCE LINES 168-172

.. code-block:: Python

    targetHSIC = ot.HSICEstimatorTargetSensitivity(
        covarianceModelCollection, X, Y, estimatorType, filterFunction
    )








.. GENERATED FROM PYTHON SOURCE LINES 173-174

We get the R2-HSIC indices:

.. GENERATED FROM PYTHON SOURCE LINES 174-178

.. code-block:: Python

    R2HSICIndices = targetHSIC.getR2HSICIndices()
    print("\n Target HSIC analysis")
    print("R2-HSIC Indices: ", R2HSICIndices)





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

 .. code-block:: none


     Target HSIC analysis
    R2-HSIC Indices:  [0.260374,0.00207302,-0.00658141]




.. GENERATED FROM PYTHON SOURCE LINES 179-180

and the HSIC indices:

.. GENERATED FROM PYTHON SOURCE LINES 180-183

.. code-block:: Python

    HSICIndices = targetHSIC.getHSICIndices()
    print("HSIC Indices: ", HSICIndices)





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

 .. code-block:: none

    HSIC Indices:  [0.00108746,8.96771e-06,-2.80515e-05]




.. GENERATED FROM PYTHON SOURCE LINES 184-185

The p-value by permutation.

.. GENERATED FROM PYTHON SOURCE LINES 185-188

.. code-block:: Python

    pvperm = targetHSIC.getPValuesPermutation()
    print("p-value (permutation): ", pvperm)





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

 .. code-block:: none

    p-value (permutation):  [0,0.237624,0.693069]




.. GENERATED FROM PYTHON SOURCE LINES 189-190

We have an asymptotic estimate of the value for this estimator.

.. GENERATED FROM PYTHON SOURCE LINES 190-193

.. code-block:: Python

    pvas = targetHSIC.getPValuesAsymptotic()
    print("p-value (asymptotic): ", pvas)





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

 .. code-block:: none

    p-value (asymptotic):  [1.42697e-09,0.316344,0.59066]




.. GENERATED FROM PYTHON SOURCE LINES 194-195

We vizualise the results.

.. GENERATED FROM PYTHON SOURCE LINES 195-208

.. code-block:: Python

    graph5 = targetHSIC.drawHSICIndices()
    view5 = otv.View(graph5)

    graph6 = targetHSIC.drawPValuesAsymptotic()
    view6 = otv.View(graph6)

    graph7 = targetHSIC.drawR2HSICIndices()
    view7 = otv.View(graph7)

    graph8 = targetHSIC.drawPValuesPermutation()
    view8 = otv.View(graph8)





.. rst-class:: sphx-glr-horizontal


    *

      .. image-sg:: /auto_reliability_sensitivity/sensitivity_analysis/images/sphx_glr_plot_hsic_estimators_ishigami_005.png
         :alt: HSIC indices
         :srcset: /auto_reliability_sensitivity/sensitivity_analysis/images/sphx_glr_plot_hsic_estimators_ishigami_005.png
         :class: sphx-glr-multi-img

    *

      .. image-sg:: /auto_reliability_sensitivity/sensitivity_analysis/images/sphx_glr_plot_hsic_estimators_ishigami_006.png
         :alt: Asymptotic p-values
         :srcset: /auto_reliability_sensitivity/sensitivity_analysis/images/sphx_glr_plot_hsic_estimators_ishigami_006.png
         :class: sphx-glr-multi-img

    *

      .. image-sg:: /auto_reliability_sensitivity/sensitivity_analysis/images/sphx_glr_plot_hsic_estimators_ishigami_007.png
         :alt: R2-HSIC indices
         :srcset: /auto_reliability_sensitivity/sensitivity_analysis/images/sphx_glr_plot_hsic_estimators_ishigami_007.png
         :class: sphx-glr-multi-img

    *

      .. image-sg:: /auto_reliability_sensitivity/sensitivity_analysis/images/sphx_glr_plot_hsic_estimators_ishigami_008.png
         :alt: p-values by permutation
         :srcset: /auto_reliability_sensitivity/sensitivity_analysis/images/sphx_glr_plot_hsic_estimators_ishigami_008.png
         :class: sphx-glr-multi-img





.. GENERATED FROM PYTHON SOURCE LINES 209-215

The Conditional HSIC estimator
------------------------------

In this last section we preprocess the input variables: that is the conditional
analysis. To do so, one has to work with a weight function.
Here the weight function is the filter function we used previously.

.. GENERATED FROM PYTHON SOURCE LINES 215-217

.. code-block:: Python

    weightFunction = filterFunction








.. GENERATED FROM PYTHON SOURCE LINES 218-219

We have to select a biased -but asymptotically unbiased- estimator

.. GENERATED FROM PYTHON SOURCE LINES 219-222

.. code-block:: Python

    estimatorType = ot.HSICVStat()









.. GENERATED FROM PYTHON SOURCE LINES 223-224

We build the conditional HSIC estimator

.. GENERATED FROM PYTHON SOURCE LINES 224-228

.. code-block:: Python

    condHSIC = ot.HSICEstimatorConditionalSensitivity(
        covarianceModelCollection, X, Y, weightFunction
    )








.. GENERATED FROM PYTHON SOURCE LINES 229-230

We get the R2-HSIC indices:

.. GENERATED FROM PYTHON SOURCE LINES 230-234

.. code-block:: Python

    R2HSICIndices = condHSIC.getR2HSICIndices()
    print("\n Conditional HSIC analysis")
    print("R2-HSIC Indices: ", R2HSICIndices)





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

 .. code-block:: none


     Conditional HSIC analysis
    R2-HSIC Indices:  [0.155438,0.017438,0.188395]




.. GENERATED FROM PYTHON SOURCE LINES 235-236

and the HSIC indices:

.. GENERATED FROM PYTHON SOURCE LINES 236-239

.. code-block:: Python

    HSICIndices = condHSIC.getHSICIndices()
    print("HSIC Indices: ", HSICIndices)





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

 .. code-block:: none

    HSIC Indices:  [0.00492342,0.000779929,0.00858311]




.. GENERATED FROM PYTHON SOURCE LINES 240-241

For the conditional estimator we only have access to the p-value by permutation:

.. GENERATED FROM PYTHON SOURCE LINES 241-245

.. code-block:: Python

    pvperm = condHSIC.getPValuesPermutation()
    print("p-value (permutation): ", pvperm)
    print("")





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

 .. code-block:: none

    p-value (permutation):  [0.029703,0.712871,0]





.. GENERATED FROM PYTHON SOURCE LINES 246-247

We can vizualise the results thanks to the various drawing methods.

.. GENERATED FROM PYTHON SOURCE LINES 247-257

.. code-block:: Python

    graph9 = condHSIC.drawHSICIndices()
    view9 = otv.View(graph9)

    graph10 = condHSIC.drawR2HSICIndices()
    view10 = otv.View(graph10)

    graph11 = condHSIC.drawPValuesPermutation()
    view11 = otv.View(graph11)





.. rst-class:: sphx-glr-horizontal


    *

      .. image-sg:: /auto_reliability_sensitivity/sensitivity_analysis/images/sphx_glr_plot_hsic_estimators_ishigami_009.png
         :alt: HSIC indices
         :srcset: /auto_reliability_sensitivity/sensitivity_analysis/images/sphx_glr_plot_hsic_estimators_ishigami_009.png
         :class: sphx-glr-multi-img

    *

      .. image-sg:: /auto_reliability_sensitivity/sensitivity_analysis/images/sphx_glr_plot_hsic_estimators_ishigami_010.png
         :alt: R2-HSIC indices
         :srcset: /auto_reliability_sensitivity/sensitivity_analysis/images/sphx_glr_plot_hsic_estimators_ishigami_010.png
         :class: sphx-glr-multi-img

    *

      .. image-sg:: /auto_reliability_sensitivity/sensitivity_analysis/images/sphx_glr_plot_hsic_estimators_ishigami_011.png
         :alt: p-values by permutation
         :srcset: /auto_reliability_sensitivity/sensitivity_analysis/images/sphx_glr_plot_hsic_estimators_ishigami_011.png
         :class: sphx-glr-multi-img





.. GENERATED FROM PYTHON SOURCE LINES 258-259

.. code-block:: Python

    otv.View.ShowAll()








.. _sphx_glr_download_auto_reliability_sensitivity_sensitivity_analysis_plot_hsic_estimators_ishigami.py:

.. only:: html

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

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

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

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

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