.. 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_quantile_confidence_chemical_process.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_quantile_confidence_chemical_process.py>`
        to download the full example code.

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

.. _sphx_glr_auto_data_analysis_sample_analysis_plot_quantile_confidence_chemical_process.py:


Estimate quantile confidence intervals from chemical process data
=================================================================

.. GENERATED FROM PYTHON SOURCE LINES 7-12

In this example, we introduce use two methods to estimate confidence intervals of the
:math:`\alpha` level quantile (:math:`\alpha \in [0,1]`).
This example is adapted from [meeker2017]_ pages 76 to 81.

See  :ref:`quantile_confidence_estimation`  and :ref:`quantile_asymptotic_confidence_estimation` to get details on the signification of these confidence interval.

.. GENERATED FROM PYTHON SOURCE LINES 12-16

.. code-block:: Python

    import openturns as ot
    import openturns.experimental as otexp
    import openturns.viewer as otv








.. GENERATED FROM PYTHON SOURCE LINES 17-19

The data represents the ordered measured amount of a chemical compound in parts per million (ppm)
taken from :math:`n=100` randomly selected batches from the process (see 5.1 data p. 76).

.. GENERATED FROM PYTHON SOURCE LINES 19-123

.. code-block:: Python

    data = [
        1.49,
        1.66,
        2.05,
        2.24,
        2.29,
        2.69,
        2.77,
        2.77,
        3.10,
        3.23,
        3.28,
        3.29,
        3.31,
        3.36,
        3.84,
        4.04,
        4.09,
        4.13,
        4.14,
        4.16,
        4.57,
        4.63,
        4.83,
        5.06,
        5.17,
        5.19,
        5.89,
        5.97,
        6.28,
        6.38,
        6.51,
        6.53,
        6.54,
        6.55,
        6.83,
        7.08,
        7.28,
        7.53,
        7.54,
        7.62,
        7.81,
        7.87,
        7.94,
        8.43,
        8.70,
        8.97,
        8.98,
        9.13,
        9.14,
        9.22,
        9.24,
        9.30,
        9.44,
        9.69,
        9.86,
        9.99,
        11.28,
        11.37,
        12.03,
        12.32,
        12.93,
        13.03,
        13.09,
        13.43,
        13.58,
        13.70,
        14.17,
        14.36,
        14.96,
        15.89,
        16.57,
        16.60,
        16.85,
        17.16,
        17.17,
        17.74,
        18.04,
        18.78,
        19.84,
        20.45,
        20.89,
        22.28,
        22.48,
        23.66,
        24.33,
        24.72,
        25.46,
        25.67,
        25.77,
        26.64,
        28.28,
        28.28,
        29.07,
        29.16,
        31.14,
        31.83,
        33.24,
        37.32,
        53.43,
        58.11,
    ]
    sample = ot.Sample.BuildFromPoint(data)








.. GENERATED FROM PYTHON SOURCE LINES 124-126

In the example, we consider the quantile of level :math:`\alpha = 10\%`,
with a confidence level of :math:`\beta = 95\%` (see example 5.7 p. 85).

.. GENERATED FROM PYTHON SOURCE LINES 126-130

.. code-block:: Python

    alpha = 0.1
    beta = 0.95
    algo = otexp.QuantileConfidence(alpha, beta)








.. GENERATED FROM PYTHON SOURCE LINES 131-142

Estimate bilateral rank: math:`(\ell,u)` such that :math:`0 \leq \ell \leq u \leq \sampleSize -1`
and defined by:

.. math::

   \begin{array}{ll}
   (\ell,u) = & \argmin \Prob{X_{(\ell)} \leq x_{\alpha} \leq X_{(u)}}\\
          & \mbox{s.t.} \Prob{X_{(\ell)} \leq x_{\alpha} \leq X_{(u)}} \geq \beta
   \end{array}

Care: indices are given in the :math:`\llbracket 0, n-1 \rrbracket` integer interval whereas the book gives them in :math:`\llbracket 1, n \rrbracket`.

.. GENERATED FROM PYTHON SOURCE LINES 142-146

.. code-block:: Python

    n = len(sample)
    l, u = algo.computeBilateralRank(n)
    print(f"l={l} u={u}")





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

 .. code-block:: none

    l=3 u=15




.. GENERATED FROM PYTHON SOURCE LINES 147-148

We can directly estimate the bilateral confidence interval of the 0.1 quantile defined by the order statistics :math:`X_{(l)}` and :math:`X_{(u)}`.

.. GENERATED FROM PYTHON SOURCE LINES 148-151

.. code-block:: Python

    ci = algo.computeBilateralConfidenceInterval(sample)
    print(f"ci={ci}")





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

 .. code-block:: none

    ci=[2.24, 4.04]




.. GENERATED FROM PYTHON SOURCE LINES 152-154

In this example, we consider the quantile of level :math:`\alpha = 90\%`,
with a confidence level of :math:`\beta = 95\%` (see example 5.1 p. 81).

.. GENERATED FROM PYTHON SOURCE LINES 154-157

.. code-block:: Python

    new_alpha = 0.90
    algo.setAlpha(new_alpha)








.. GENERATED FROM PYTHON SOURCE LINES 158-159

We get the bilateral confidence interval of the 0.91 quantile.

.. GENERATED FROM PYTHON SOURCE LINES 159-162

.. code-block:: Python

    l, u = algo.computeBilateralRank(n)
    print(f"l={l} u={u}")





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

 .. code-block:: none

    l=84 u=96




.. GENERATED FROM PYTHON SOURCE LINES 163-166

.. code-block:: Python

    ci = algo.computeBilateralConfidenceInterval(sample)
    print(f"ci={ci}")





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

 .. code-block:: none

    ci=[24.33, 33.24]




.. GENERATED FROM PYTHON SOURCE LINES 167-175

We can estimate the rank :math:`k_{low}` which is the largest rank :math:`k`  with
:math:`0 \leq k \leq \sampleSize -1` such that:

.. math::

   \Prob{X_{(k)} \leq x_{\alpha}} \geq \beta.

We notice that the order statistics of the lower bound is the same as in the bilateral confidence interval.

.. GENERATED FROM PYTHON SOURCE LINES 175-178

.. code-block:: Python

    k_low = algo.computeUnilateralRank(n, True)
    print(f"k_low={k_low}")





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

 .. code-block:: none

    k_low=84




.. GENERATED FROM PYTHON SOURCE LINES 179-182

In other words, the interval :math:`\left[ X_{(k_{low})}, +\infty\right[` is a unilateral
confidence interval for the 0.9 quantile with confidence :math:`\beta`. We can directly
estimate this interval.

.. GENERATED FROM PYTHON SOURCE LINES 182-185

.. code-block:: Python

    ci_low = algo.computeUnilateralConfidenceInterval(sample, True)
    print(f"ci_low={ci_low}")





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

 .. code-block:: none

    ci_low=[24.33, (1.79769e+308) +inf[




.. GENERATED FROM PYTHON SOURCE LINES 186-194

We now estimate the rank :math:`k_{up}` which is the smallest rank :math:`k`  with
:math:`0 \leq k \leq \sampleSize -1` such that:

.. math::

   \Prob{X_{(k)} \geq x_{\alpha}} \geq \beta.

We notice that the order statistics of the upper bound is slightly smaller than in the bilateral confidence interval.

.. GENERATED FROM PYTHON SOURCE LINES 194-197

.. code-block:: Python

    k_up = algo.computeUnilateralRank(n)
    print(f"k_up={k_up}")





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

 .. code-block:: none

    k_up=95




.. GENERATED FROM PYTHON SOURCE LINES 198-201

In other words, the interval :math:`\left]-\infty,  X_{(k_{up})}\right]` is a unilateral
confidence interval for the 0.9 quantile with confidence :math:`\beta`. We can directly
estimate this interval.

.. GENERATED FROM PYTHON SOURCE LINES 201-204

.. code-block:: Python

    ci_up = algo.computeUnilateralConfidenceInterval(sample)
    print(f"ci_up={ci_up}")





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

 .. code-block:: none

    ci_up=]-inf (-1.79769e+308), 31.83]




.. GENERATED FROM PYTHON SOURCE LINES 205-206

We had the empirical estimation of the quantile, with the order statistics :math:`X_{(\lfloor \sampleSize \alpha, \rfloor)}`.

.. GENERATED FROM PYTHON SOURCE LINES 206-208

.. code-block:: Python

    emp_quant = sample.computeQuantile(new_alpha)[0]








.. GENERATED FROM PYTHON SOURCE LINES 209-212

We illustrate here the confidence intervals we obtained. To do that, we
draw the empirical cumulative distribution function and the bounds of the bilateral confidence intervals.
We first draw the empirical cumulative distribution function.

.. GENERATED FROM PYTHON SOURCE LINES 212-216

.. code-block:: Python

    user_defined_dist = ot.UserDefined(sample)
    g = user_defined_dist.drawCDF(sample.getMin(), sample.getMax() * 1.5)
    g = user_defined_dist.drawCDF(0.0, 60.0)








.. GENERATED FROM PYTHON SOURCE LINES 217-219

Then we had the bounds of the bilateral confidence intervals.
First the bilateral interval.

.. GENERATED FROM PYTHON SOURCE LINES 219-225

.. code-block:: Python

    line_bil_low = ot.Curve(
        [ci.getLowerBound(), ci.getLowerBound()],
        [[-0.15], [user_defined_dist.computeCDF(ci.getLowerBound())]],
    )
    line_bil_low.setLineStyle("dashed")








.. GENERATED FROM PYTHON SOURCE LINES 226-233

.. code-block:: Python

    line_bil_up = ot.Curve(
        [ci.getUpperBound(), ci.getUpperBound()],
        [[-0.20], [user_defined_dist.computeCDF(ci.getUpperBound())]],
    )
    line_bil_up.setLineStyle("dashed")
    line_bil_up.setColor(line_bil_low.getColor())








.. GENERATED FROM PYTHON SOURCE LINES 234-237

.. code-block:: Python

    g.add(line_bil_low)
    g.add(line_bil_up)








.. GENERATED FROM PYTHON SOURCE LINES 238-239

Then the unilateral confidence intervals.

.. GENERATED FROM PYTHON SOURCE LINES 239-245

.. code-block:: Python

    line_unil_low = ot.Curve(
        [ci_low.getLowerBound(), ci_low.getLowerBound()],
        [[0.0], [user_defined_dist.computeCDF(ci_low.getLowerBound())]],
    )
    line_unil_low.setLineStyle("dotted")








.. GENERATED FROM PYTHON SOURCE LINES 246-252

.. code-block:: Python

    line_unil_up = ot.Curve(
        [ci_up.getUpperBound(), ci_up.getUpperBound()],
        [[-0.05], [user_defined_dist.computeCDF(ci_up.getUpperBound())]],
    )
    line_unil_up.setLineStyle("dashed")








.. GENERATED FROM PYTHON SOURCE LINES 253-256

.. code-block:: Python

    g.add(line_unil_low)
    g.add(line_unil_up)








.. GENERATED FROM PYTHON SOURCE LINES 257-258

At last, the empirical estimation.

.. GENERATED FROM PYTHON SOURCE LINES 258-268

.. code-block:: Python

    line_emp = ot.Curve(
        [[emp_quant], [emp_quant], [0.0]],
        [
            [-0.10],
            [user_defined_dist.computeCDF(emp_quant)],
            [user_defined_dist.computeCDF(emp_quant)],
        ],
    )
    line_bil_up.setLineStyle("dashed")








.. GENERATED FROM PYTHON SOURCE LINES 269-271

.. code-block:: Python

    g.add(line_emp)








.. GENERATED FROM PYTHON SOURCE LINES 272-273

We add some labels.

.. GENERATED FROM PYTHON SOURCE LINES 273-286

.. code-block:: Python

    text_bil_low = ot.Text([[ci.getLowerBound()[0], -0.2]], ["bilat. low. bound"])
    text_bil_up = ot.Text([[ci.getUpperBound()[0], -0.25]], ["bilat. up. bound"])
    text_low = ot.Text([[ci_low.getLowerBound()[0], -0.05]], ["unilat. low. bound"])
    text_up = ot.Text([[ci_up.getUpperBound()[0], -0.1]], ["unilat. up. bound"])
    text_emp = ot.Text([[emp_quant, -0.15]], ["emp quant"])
    text_emp_2 = ot.Text([[-0.1, 0.9]], ["quant level"])
    g.add(text_bil_low)
    g.add(text_bil_up)
    g.add(text_low)
    g.add(text_up)
    g.add(text_emp)
    g.add(text_emp_2)








.. GENERATED FROM PYTHON SOURCE LINES 287-320

.. code-block:: Python

    g.setColors(
        [
            "#1f77b4",
            "#ff7f0e",
            "#ff7f0e",
            "#9467bd",
            "#d62728",
            "black",
            "#ff7f0e",
            "#ff7f0e",
            "#9467bd",
            "#d62728",
            "black",
            "black",
        ]
    )
    g.setLegends(
        [
            "Emp. CDF",
            "Bilat CI: lower bound",
            "Bilat CI: upper bound",
            "Unilat CI: lower bound",
            "Unilat CI: upper bound",
            "Emp quant",
        ]
    )

    g.setLegendCorner([1.0, 1.0])
    g.setLegendPosition("upper left")
    g.setTitle("Estimation of the quantile of level 0.9")
    g.setXTitle("x")
    view = otv.View(g)




.. image-sg:: /auto_data_analysis/sample_analysis/images/sphx_glr_plot_quantile_confidence_chemical_process_001.svg
   :alt: Estimation of the quantile of level 0.9
   :srcset: /auto_data_analysis/sample_analysis/images/sphx_glr_plot_quantile_confidence_chemical_process_001.svg
   :class: sphx-glr-single-img





.. GENERATED FROM PYTHON SOURCE LINES 321-322

Display all the figures.

.. GENERATED FROM PYTHON SOURCE LINES 322-323

.. code-block:: Python

    otv.View.ShowAll()








.. _sphx_glr_download_auto_data_analysis_sample_analysis_plot_quantile_confidence_chemical_process.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_confidence_chemical_process.ipynb <plot_quantile_confidence_chemical_process.ipynb>`

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

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

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

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