.. 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_estimation.py" .. LINE NUMBERS ARE GIVEN BELOW. .. only:: html .. note:: :class: sphx-glr-download-link-note :ref:`Go to the end ` to download the full example code. .. rst-class:: sphx-glr-example-title .. _sphx_glr_auto_data_analysis_sample_analysis_plot_quantile_confidence_estimation.py: Estimate quantile confidence intervals ====================================== .. GENERATED FROM PYTHON SOURCE LINES 7-19 In this example, we introduce two methods to estimate a confidence interval of the :math:`\alpha` level quantile (:math:`\alpha \in [0,1]`) of the distribution of a random variable :math:`X` of dimension 1. Both methods use the order statistics to estimate: - an asymptotic confidence interval with confidence level :math:`\beta \in [0,1]`, - an exact upper bounded confidence interval with confidence level :math:`\beta \in [0,1]`. In this example, we consider the quantile of level :math:`\alpha = 95\%`, with a confidence level of :math:`\beta = 90\%`. 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 19-24 .. code-block:: Python import openturns as ot import openturns.experimental as otexp import math as m .. GENERATED FROM PYTHON SOURCE LINES 25-26 We consider a random vector which is the output of a model and an input distribution. .. GENERATED FROM PYTHON SOURCE LINES 26-33 .. 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) .. GENERATED FROM PYTHON SOURCE LINES 34-35 Create the output random vector .. GENERATED FROM PYTHON SOURCE LINES 35-37 .. code-block:: Python output = ot.CompositeRandomVector(model, inputVector) .. GENERATED FROM PYTHON SOURCE LINES 38-39 We define the level :math:`\alpha` of the quantile and the confidence level :math:`\beta`. .. GENERATED FROM PYTHON SOURCE LINES 39-42 .. code-block:: Python alpha = 0.95 beta = 0.90 .. GENERATED FROM PYTHON SOURCE LINES 43-44 We generate a sample of the variable. .. GENERATED FROM PYTHON SOURCE LINES 44-47 .. code-block:: Python n = 10**4 sample = output.getSample(n) .. GENERATED FROM PYTHON SOURCE LINES 48-51 We get the empirical estimator of the :math:`\alpha` level quantile which is the :math:`\lfloor \sampleSize \alpha \rfloor` -th order statistics evaluated on the sample. .. GENERATED FROM PYTHON SOURCE LINES 51-54 .. code-block:: Python empiricalQuantile = sample.computeQuantile(alpha) print(empiricalQuantile) .. rst-class:: sphx-glr-script-out .. code-block:: none [4.31123] .. GENERATED FROM PYTHON SOURCE LINES 55-71 The asymptotic confidence interval of level :math:`\beta` is :math:`\left[ X_{(i_n)}, X_{(j_n)}\right]` such that: .. math:: i_\sampleSize & = \left\lfloor \sampleSize \alpha - \sqrt{\sampleSize} \; z_{\frac{1+\beta}{2}} \; \sqrt{\alpha(1 - \alpha)} \right\rfloor\\ j_\sampleSize & = \left\lfloor \sampleSize \alpha + \sqrt{\sampleSize} \; z_{\frac{1+\beta}{2}} \; \sqrt{\alpha(1 - \alpha)} \right\rfloor where :math:`z_{\frac{1+\beta}{2}}` is the :math:`\frac{1+\beta}{2}` level quantile of the standard normal distribution (see [delmas2006]_ proposition 11.1.13). Then we have: .. math:: \lim\limits_{\sampleSize \rightarrow +\infty} \Prob{x_{\alpha} \in \left[ X_{(i_\sampleSize,\sampleSize)}, X_{(j_\sampleSize,\sampleSize)}\right]} = \beta .. GENERATED FROM PYTHON SOURCE LINES 71-76 .. code-block:: Python a_beta = ot.Normal(1).computeQuantile((1.0 + beta) / 2.0)[0] i_n = int(n * alpha - a_beta * m.sqrt(n * alpha * (1.0 - alpha))) j_n = int(n * alpha + a_beta * m.sqrt(n * alpha * (1.0 - alpha))) print(i_n, j_n) .. rst-class:: sphx-glr-script-out .. code-block:: none 9464 9535 .. GENERATED FROM PYTHON SOURCE LINES 77-78 Get the sorted sample .. GENERATED FROM PYTHON SOURCE LINES 78-80 .. code-block:: Python sortedSample = sample.sort() .. GENERATED FROM PYTHON SOURCE LINES 81-85 Get the asymptotic confidence interval :math:`\left[ X_{(i_n)}, X_{(j_n)}\right]`. Care: the index in the sorted sample is :math:`i_n-1` and :math:`j_n-1` because python numbering starts at 0. .. GENERATED FROM PYTHON SOURCE LINES 85-89 .. code-block:: Python infQuantile = sortedSample[i_n - 1] supQuantile = sortedSample[j_n - 1] print(infQuantile, empiricalQuantile, supQuantile) .. rst-class:: sphx-glr-script-out .. code-block:: none [4.20242] [4.31123] [4.4304] .. GENERATED FROM PYTHON SOURCE LINES 90-92 This can be done using :class:`~openturns.experimental.QuantileConfidence` when the ranks :math:`i_n` and :math:`j_n` are directly given in :math:`\llbracket 0, \sampleSize-1 \rrbracket`. .. GENERATED FROM PYTHON SOURCE LINES 92-97 .. code-block:: Python algo = otexp.QuantileConfidence(alpha, beta) i_n, j_n = algo.computeAsymptoticBilateralRank(n) ci = algo.computeAsymptoticBilateralConfidenceInterval(sample) print(f"asymptotic ranks={[i_n, j_n]} CI={ci}") .. rst-class:: sphx-glr-script-out .. code-block:: none asymptotic ranks=[9463, 9534] CI=[4.20242, 4.4304] .. GENERATED FROM PYTHON SOURCE LINES 98-102 The empirical quantile was estimated with the :math:`\lfloor \sampleSize\alpha \rfloor` -th order statistics evaluated on the sample of size :math:`\sampleSize`. We define :math:`i = \sampleSize-\lfloor \sampleSize\alpha \rfloor` and we evaluate the minimum sample size :math:`\tilde{\sampleSize}` that ensures that the :math:`(\tilde{\sampleSize}-i)` order statistics is greater than :math:`x_{\alpha}` with the confidence :math:`\beta`. .. GENERATED FROM PYTHON SOURCE LINES 102-106 .. code-block:: Python i = n - int(n * alpha) minSampleSize = algo.computeUnilateralMinimumSampleSize(i, True) print(minSampleSize) .. rst-class:: sphx-glr-script-out .. code-block:: none 10583 .. GENERATED FROM PYTHON SOURCE LINES 107-109 Now we can evaluate the upper bounded confidence interval: we generate a sample with the previous minimum size and extract the empirical quantile of order :math:`(\tilde{\sampleSize}-i)`. .. GENERATED FROM PYTHON SOURCE LINES 109-113 .. code-block:: Python sample = output.getSample(minSampleSize) upperBoundQuantile = sample.sort()[-i - 1] print(upperBoundQuantile) .. rst-class:: sphx-glr-script-out .. code-block:: none [4.44384] .. GENERATED FROM PYTHON SOURCE LINES 114-115 We can also find this rank back from the size. .. GENERATED FROM PYTHON SOURCE LINES 115-118 .. code-block:: Python k = algo.computeUnilateralRank(minSampleSize, True) print(k, minSampleSize - i - 1) .. rst-class:: sphx-glr-script-out .. code-block:: none 10024 10082 .. GENERATED FROM PYTHON SOURCE LINES 119-120 The quantile bound is given in the same manner from the sample and the rank. .. GENERATED FROM PYTHON SOURCE LINES 120-122 .. code-block:: Python ci = algo.computeUnilateralConfidenceInterval(sample, True) print(ci) .. rst-class:: sphx-glr-script-out .. code-block:: none [4.25546, (1.79769e+308) +inf[ .. _sphx_glr_download_auto_data_analysis_sample_analysis_plot_quantile_confidence_estimation.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_estimation.ipynb ` .. container:: sphx-glr-download sphx-glr-download-python :download:`Download Python source code: plot_quantile_confidence_estimation.py ` .. container:: sphx-glr-download sphx-glr-download-zip :download:`Download zipped: plot_quantile_confidence_estimation.zip `