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

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

.. _sphx_glr_auto_data_analysis_manage_data_and_samples_plot_quantile_estimation_wilks.py:


Estimate Wilks and empirical quantile
=====================================

.. GENERATED FROM PYTHON SOURCE LINES 7-47

In this example we want to evaluate a particular quantile, with the empirical estimator or the Wilks one, from a sample of a random variable.


Let us suppose we want to estimate the quantile :math:`q_{\alpha}` of order :math:`\alpha` of the variable :math:`Y`:
:math:`P(Y \leq q_{\alpha}) = \alpha`, from the sample :math:`(Y_1, ..., Y_n)`
of size :math:`n`, with a confidence level equal to :math:`\beta`.

We note :math:`(Y^{(1)}, ..., Y^{(n)})` the sample where the values are sorted in ascending order.
The empirical estimator, noted :math:`q_{\alpha}^{emp}`, and its confidence interval, is defined by the expressions:

.. math::
    \left\{
     \begin{array}{lcl}
       q_{\alpha}^{emp} & = & Y^{(E[n\alpha])} \\
       P(q_{\alpha} \in [Y^{(i_n)}, Y^{(j_n)}]) & = & \beta \\
       i_n & = & E[n\alpha - a_{\alpha}\sqrt{n\alpha(1-\alpha)}] \\
       i_n & = & E[n\alpha + a_{\alpha}\sqrt{n\alpha(1-\alpha)}]
     \end{array}
    \right\}

The Wilks estimator, noted :math:`q_{\alpha, \beta}^{Wilks}`, and its confidence interval, is defined by the expressions:

.. math::
  \left\{
  \begin{array}{lcl}
    q_{\alpha, \beta}^{Wilks} & = & Y^{(n-i)} \\
    P(q_{\alpha}  \leq q_{\alpha, \beta}^{Wilks}) & \geq & \beta \\
    i\geq 0 \, \, /  \, \, n \geq N_{Wilks}(\alpha, \beta,i)
  \end{array}
  \right\}

Once the order :math:`i` has been chosen, the Wilks number :math:`N_{Wilks}(\alpha, \beta,i)` is evaluated,
thanks to the static method :math:`ComputeSampleSize(\alpha, \beta, i)` of the Wilks object.

In the example, we want to evaluate a quantile :math:`\alpha = 95\%`,
with a confidence level of :math:`\beta = 90\%` thanks to the :math:`4`th maximum of
the ordered sample (associated to the order :math:`i = 3`).

Be careful: :math:`i=0` means that the Wilks estimator is the maximum of the sample:
it corresponds to the first maximum of the sample.

.. GENERATED FROM PYTHON SOURCE LINES 49-55

.. code-block:: Python

    import openturns as ot
    import math as m
    import openturns.viewer as viewer

    ot.Log.Show(ot.Log.NONE)








.. GENERATED FROM PYTHON SOURCE LINES 56-66

.. code-block:: Python

    model = ot.SymbolicFunction(["x1", "x2"], ["x1^2+x2"])
    R = ot.CorrelationMatrix(2)
    R[0, 1] = -0.6
    inputDist = ot.Normal([0.0, 0.0], R)
    inputDist.setDescription(["X1", "X2"])
    inputVector = ot.RandomVector(inputDist)

    # Create the output random vector Y=model(X)
    output = ot.CompositeRandomVector(model, inputVector)








.. GENERATED FROM PYTHON SOURCE LINES 67-68

Quantile level

.. GENERATED FROM PYTHON SOURCE LINES 68-73

.. code-block:: Python

    alpha = 0.95

    # Confidence level of the estimation
    beta = 0.90








.. GENERATED FROM PYTHON SOURCE LINES 74-75

Get a sample of the variable

.. GENERATED FROM PYTHON SOURCE LINES 75-80

.. code-block:: Python

    N = 10**4
    sample = output.getSample(N)
    graph = ot.UserDefined(sample).drawCDF()
    view = viewer.View(graph)




.. image-sg:: /auto_data_analysis/manage_data_and_samples/images/sphx_glr_plot_quantile_estimation_wilks_001.png
   :alt: y0 CDF
   :srcset: /auto_data_analysis/manage_data_and_samples/images/sphx_glr_plot_quantile_estimation_wilks_001.png
   :class: sphx-glr-single-img





.. GENERATED FROM PYTHON SOURCE LINES 81-82

Empirical Quantile Estimator

.. GENERATED FROM PYTHON SOURCE LINES 82-98

.. code-block:: Python

    empiricalQuantile = sample.computeQuantile(alpha)

    # Get the indices of the confidence interval bounds
    aAlpha = ot.Normal(1).computeQuantile((1.0 + beta) / 2.0)[0]
    min_i = int(N * alpha - aAlpha * m.sqrt(N * alpha * (1.0 - alpha)))
    max_i = int(N * alpha + aAlpha * m.sqrt(N * alpha * (1.0 - alpha)))
    # print(min_i, max_i)

    # Get the sorted sample
    sortedSample = sample.sort()

    # Get the Confidence interval of the Empirical Quantile Estimator [infQuantile, supQuantile]
    infQuantile = sortedSample[min_i - 1]
    supQuantile = sortedSample[max_i - 1]
    print(infQuantile, empiricalQuantile, supQuantile)





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

 .. code-block:: none

    [4.13903] [4.28037] [4.35925]




.. GENERATED FROM PYTHON SOURCE LINES 99-100

Wilks number

.. GENERATED FROM PYTHON SOURCE LINES 100-104

.. code-block:: Python

    i = N - (min_i + max_i) // 2  # compute wilks with the same sample size
    wilksNumber = ot.Wilks.ComputeSampleSize(alpha, beta, i)
    print("wilksNumber =", wilksNumber)





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

 .. code-block:: none

    wilksNumber = 10604




.. GENERATED FROM PYTHON SOURCE LINES 105-106

Wilks Quantile Estimator

.. GENERATED FROM PYTHON SOURCE LINES 106-109

.. code-block:: Python

    algo = ot.Wilks(output)
    wilksQuantile = algo.computeQuantileBound(alpha, beta, i)
    print("wilks Quantile 0.95 =", wilksQuantile)




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

 .. code-block:: none

    wilks Quantile 0.95 = [4.37503]





.. _sphx_glr_download_auto_data_analysis_manage_data_and_samples_plot_quantile_estimation_wilks.py:

.. only:: html

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

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

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

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

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