scipy.stats._quantile:quantile

Notes

Given a sample x from an underlying distribution, quantile provides a nonparametric estimate of the inverse cumulative distribution function.

By default, this is done by interpolating between adjacent elements in y, a sorted copy of x:

(1-g)*y[j] + g*y[j+1]

where the index j and coefficient g are the integral and fractional components of p * (n-1), and n is the number of elements in the sample.

This is a special case of Equation 1 of H&F ^[1]. More generally,

• j = (p*n + m - 1) // 1, and

• g = (p*n + m - 1) % 1,

where m may be defined according to several different conventions. The preferred convention may be selected using the method parameter:

=============================== =============== ===============
``method``                      number in H&F   ``m``
=============================== =============== ===============
``interpolated_inverted_cdf``   4               ``0``
``hazen``                       5               ``1/2``
``weibull``                     6               ``p``
``linear`` (default)            7               ``1 - p``
``median_unbiased``             8               ``p/3 + 1/3``
``normal_unbiased``             9               ``p/4 + 3/8``
=============================== =============== ===============

Note that indices j and j + 1 are clipped to the range 0 to n - 1 when the results of the formula would be outside the allowed range of non-negative indices. When j is clipped to zero, g is set to zero as well. The -1 in the formulas for j and g accounts for Python's 0-based indexing.

The table above includes only the estimators from ^[1] that are continuous functions of probability p (estimators 4-9). SciPy also provides the three discontinuous estimators from ^[1] (estimators 1-3), where j is defined as above, m is defined as follows, and g is 0 when index = p*n + m - 1 is less than 0 and otherwise is defined below.

• inverted_cdf: m = 0 and g = int(index - j > 0)

• averaged_inverted_cdf: m = 0 and g = (1 + int(index - j > 0)) / 2

• closest_observation: m = -1/2 and g = 1 - int((index == j) & (j%2 == 1))

Note that for methods inverted_cdf and averaged_inverted_cdf, only the relative proportions of tied observations (and relative weights) affect the results; for all other methods, the total number of observations (and absolute weights) matter.

A different strategy for computing quantiles from ^[2], method='harrell-davis', uses a weighted combination of all elements. The weights are computed as:

where $n$ is the number of elements in the sample, $i$ are the indices $1,2,...,n-1,n$ of the sorted elements, $a=p(n+1)$, $b=(1-p)(n+1)$, $p$ is the probability of the quantile, and $I$ is the regularized, lower incomplete beta function (scipy.special.betainc).

method='round_nearest' is equivalent to indexing y[j], where

j = int(np.round(p*n) if p < 0.5 else np.round(n*p - 1))

This is useful when winsorizing data: replacing p*n of the most extreme observations with the next most extreme observation. method='round_outward' adjusts the direction of rounding to winsorize fewer elements

j = int(np.floor(p*n) if p < 0.5 else np.ceil(n*p - 1))

and method='round_inward' rounds to winsorize more elements

j = int(np.ceil(p*n) if p < 0.5 else np.floor(n*p - 1))

These methods are also useful for trimming data: removing p*n of the most extreme observations. See outliers for example applications.

Array API Standard Support

quantile has experimental support for Python Array API Standard compatible backends in addition to NumPy. Please consider testing these features by setting an environment variable SCIPY_ARRAY_API=1 and providing CuPy, PyTorch, JAX, or Dask arrays as array arguments. The following combinations of backend and device (or other capability) are supported.

====================  ====================  ====================
Library               CPU                   GPU
====================  ====================  ====================
NumPy                 ✅                     n/a                 
CuPy                  n/a                   ✅                   
PyTorch               ✅                     ✅                   
JAX                   ⚠️ no JIT             ⚠️ no JIT           
Dask                  ⛔                     n/a                 
====================  ====================  ====================

See dev-arrayapi for more information.