Numpy implementation of the Quantile function [duplicate]

Question

I am trying to wrap my head around the implementation of the quantile function for finite samples and, specifically, in numpy (main reason to do this: I am working on conformal predictions). It is likely I am missing something clear and welcome anyone who can point out my clear mistake!

Assume we observe values $\mathcal{X} = \{x_1, x_2 \ldots x_n\}$ and we want to calculate the $p \in [0, 1]$ quantile from the observed sample. So, we denote the quantile function over this finite sample as $Q(q, \mathcal{X})$.

If I check the Wikipedia page of quantile (https://en.wikipedia.org/wiki/Quantile), we have that the quantile are defined as follows:

$Q(p, \mathcal{X}) := x: P(X<x) \leq p \text{ and } P(X\leq x) \geq p$

I am not 100% sure how to interpret "$P(X<x)$" here. I assume this is a reference to the e.c.d.f. (https://en.wikipedia.org/wiki/Empirical_distribution_function) of $\mathcal{X}$. That is, "$P(X<x)$" is operationalized as the proportion of points in the sample less than $x$. Therefore, we get the following equivalent definition:

$Q(p, \mathcal{X}) = \text{sup}\{x \in \mathbb{R} : \text{ecdf}_\mathcal{X}(x) < p \}$.

So far so good, I believe. Assume $\mathcal{X} = \{2,3,4,10 \}$. Now let's check numpy implementation:

import numpy as np
v = np.array([2,3,4,10])
q_third  = np.quantile(v, 1/3) # 3
assert q_third == 3 # This is true! quantile is exactly 3

I am at loss here. numpy tells me that $Q(\frac{1}{3}, \{2,3,4,10 \}) = 3$. Based on the definition above, I expected that $Q(\frac{1}{4}, \{2,3,4,10 \}) = 3$ . In fact, the definition above is ambiguous w.r.t. to $ Q(\frac{1}{3}, \{2,3,4,10 \})$ as far as I understand.

Under to linear interpolation (the default method in numpy) we can assume that $Q(\frac{1}{4}, \{2,3,4,10 \}) \neq Q(\frac{1}{3}, \{2,3,4,10 \})$. So, either me or numpy are wrong. I tend to think numpy must be correct here, but I cannot spot my error. If anyone can explain to me, I would be grateful.

I acknowledge that there is a similar question (Definition of quantile), but the answer to said question simply reports the same Wikipedia definition above and refers to R (irrelevant in this context). An explanation of the logic numpy applies in my example (hopefully in the same notation I use) will be accepted as an answer.

Answer 1

The value returned by the np.quantile function will depend on the interpolation behavior; there are a number of different ways to define "quantile", and numpy implements several of them. For more detail, see R. J. Hyndman and Y. Fan, “Sample quantiles in statistical packages,” The American Statistician, 50(4), pp. 361-365, 1996.

This is not much different from the situation as arises when estimating the median for a sample with an even sample size. See also: Definition of quantile

How np.quantile does interpolation is explained in the numpy documentation.

Given a vector V of length n, the q-th quantile of V is the value q of the way from the minimum to the maximum in a sorted copy of V. The values and distances of the two nearest neighbors as well as the method parameter will determine the quantile if the normalized ranking does not match the location of q exactly. This function is the same as the median if q=0.5, the same as the minimum if q=0.0 and the same as the maximum if q=1.0.

The optional method parameter specifies the method to use when the desired quantile lies between two indexes i and j = i + 1. In that case, we first determine i + g, a virtual index that lies between i and j, where i is the floor and g is the fractional part of the index. The final result is, then, an interpolation of a[i] and a[j] based on g. During the computation of g, i and j are modified using correction constants alpha and beta whose choices depend on the method used. Finally, note that since Python uses 0-based indexing, the code subtracts another 1 from the index internally.

The following formula determines the virtual index i + g, the location of the quantile in the sorted sample: $$ i + g = q \times (n - \alpha - \beta +1 ) + \alpha $$

This is straightforward to implement; this implementation matches numpy for method="linear", which is the default.

import numpy as np


def my_quantile(a, q):
    n = a.size
    alpha = 1.0
    beta = 1.0
    i = int(q * n)
    g = q * (n - alpha - beta + 1.0) + alpha - i
    a = np.sort(a)
    upper = a[i]
    lower = a[i - 1]
    return (upper - lower) * g + lower


if __name__ == "__main__":
    v = np.array([2, 3, 4, 10])
    print(v)
    q_val = np.quantile(v, 1.0 / 3.0, method="linear")
    me = my_quantile(a=v, q=1.0 / 3.0)
    print(f"numpy:\t{q_val}") # 3.0
    print(f"mine:\t{me}") # 3.0

    q_val = np.quantile(v, 0.25, method="linear")
    me = my_quantile(a=v, q=0.25)
    print(20 * "-")
    print(f"numpy:\t{q_val}") # 2.75
    print(f"mine:\t{me}") # 2.75

    prng = np.random.default_rng(42)
    norm_v = prng.normal(size=100)
    me = my_quantile(a=norm_v, q=0.25)
    print(20 * "-")
    print(f"numpy:\t{q_val}") # -0.5439015088644885
    print(f"mine:\t{me}") # -0.5439015088644885

Uses numpy version 1.24.3.

In comments, OP writes that they expect np.quantile(v, 0.25) to return 2. This is achieved by a different method than the default ”linear” method. If OP desires this behavior, then they should use np.quantile(v, 0.25, method='lower'), which returns 2.