Let $Z_{(1)} \le \ldots \le Z_{(n)}$ denote the order statistics of a collection of independent random variables, $Z_{i}$ indexed by $i \in N$ where $|N| = n$. Let $F_i$ and $f_i$ be the distribution and density functions for each $i$.

To generate, the order statistics we do the following: $$ \begin{align*} & Z_{i} \sim F_i\quad\textbf{for } i \in N \\ & Z_{(1)}, \ldots, Z_{(n)} = \mathrm{sort}\left([ Z_{i} \textbf{ for }i \in N ] \right) \end{align*} $$

In this post, we work out the distribution and density functions for each $Z_{(k)}$.

There is no shortage of lecture notes and books online that describe the distribution functions for $Z_{(1)}$ (min) and $Z_{(n)}$ (max) as these cases are particularly simple and more commonly used; similarly, for the case of identically distributed random variables (where the $F_i$s are equal). The more general (non-identical, but still independent) case turns out to not be too much harder. Why do I care about this problem? These distribution functions arise when analyzing order-based sampling without replacement designs, such as the Gumbel-top-k design that I have written about many times.

Before continuing, I should note that my derivation and implementation draws heavily on Aires (1999); specifically, Lemma 2.

Without further ado...

Distribution Function¶

$$ \begin{align*} F^N_{(k)}(t) &\overset{\mathrm{def}}{=} \mathrm{Pr}\left( Z_{(k)} \le t \right) \\ &= \mathrm{Pr}\left( \text{at least } k \text{ of } Z_1 \ldots Z_n \text{ are } \le t \right) \\ &= \sum_{\substack{ Y \subseteq N \\ |Y| \ge k}} \underbrace{\left( \prod_{i \in Y} F_i(t) \right)}_{ \text{ values } \le t} \, \underbrace{\left( \prod_{i \in N \smallsetminus Y } \left( 1-F_i(t) \right) \right)}_{ \text{ values } > t} \end{align*} $$ It is a sum over all the ways in which $t$ can be the $k^{\text{th}}$ order statistics: there must be at least $k$ values that are $\le t$, and the remaining values that are $> t$. These events are mutually exclusive and exhaustive. Let $Y$ be the set of variables that are $\le t$. Clearly, we require $|Y| \ge k$. Because the $Z_i$s are independent, the probability of the event that $Y$ is the set of $k$ smallest elements is simply the product of the $F_i$s in $Y$ and $(1-F_i)$s in the complement of $Y$.

It turns out that we can evaluate the exponential-size sum in $\mathcal{O}(n^2)$ time with a simple algorithm, which is essentially a probabilistic generalization of the dynamic program for evaluating the binomial coefficients, ${n \choose k}$. Furthermore, we can simultaneously evaluate the distribution function for all $n$ order statistics in $\mathcal{O}(n^2)$ time. The algorithm for doing this computation efficiently is in the code section of this article.

Density Function¶

\begin{align} f^N_{(k)}(t) &\overset{\mathrm{def}}{=} \frac{\partial}{\partial t} F^N_{(k)}(t) \\ &= \frac{\partial}{\partial t} \left[ \sum_{\substack{ Y \subseteq N: \\ |Y| \ge k}} \left( \prod_{i \in Y} F_i(t) \right) \left( \prod_{i \in N \smallsetminus Y } \left(1-F_i(t)\right) \right) \right] \\ &= \sum_j \sum_{\substack{ Y \subseteq N: \\ |Y| \ge k}} (-1)^{1[j \notin Y]} f_j(t) \left( \prod_{\substack{i \in Y \\ i \ne j}} F_i(t) \right) \left( \prod_{\substack{i \in N \smallsetminus Y \\ i \ne j}} \left(1-F_i(t)\right) \right) \\ &= \sum_j \left(F^{N - j}_{(k-1)}(t) - F^{N - j}_{(k)}(t)\right) f_j(t) \end{align} where $F^{N - j}_{(k-1)}$ and $F^{N - j}_{(k)}$ are the distribution functions for the order statistics $k-1$ and $k$ over the set of elements excluding the random variable $Z_j$.

Code¶

Utility Functions to Efficiently Sum over Weighted Subsets.¶

Definition: Let $p_i$ and $q_i$ be inclusion and exclusion weights (respectively) for some set of elements $i \in N$. The total weight of subsets of $N$ of size with size least $K$ is $$ C^N_K \overset{\mathrm{def}}{=} \sum_{\substack{Y \subseteq N:\\ |Y|\ge K}} \left(\prod_{i \in Y} p_i \right)\left(\prod_{i \in N \smallsetminus Y} q_i\right) $$ We can evaluate $C^N_K$ efficiently in $\mathcal{O}(n^2)$ with the following algorithms.

The method below is used to efficiently evaluate the gradient of the C function in time $\mathcal{O}(n^2)$. This method is used to evaluated the PDF. It was derived by manually applying algorithmic differentiation, but it could also be evaluated using an automatic differentiation toolkit (read more).

Tests¶

Synethic Data

Below, we test out our analytical method for computing the distribution function of independent normal distributions. The method should work for any collection of distributions in scipy.stats. You can even mix and match the distribution types. We will compare the analytical $F_{(k)}$ to the empirical distribution of $Z_{(k)}$ using a large sample size.

CDF

Below, I use my arsenal library to fit the empirical distribution function and compare it to the analytical solution.

In the plots above, the red line is the empirical distribution funciton, and the blue dotted line is the analytical distribution function. In all cases, we see that the fit is good.

PDF

Below, we check the PDF against histograms.

Closing¶

I would like to thank Lin Zheng and the commenter "Mechanical Turk" for pointing out a mistake in the original version of this article. I would also like to thank Ryan Cotterell for valuable feedback and discussions that greatly improved this article.

This Jupyter notebook is available for download.

I am very happy to answer questions! There are many ways to do that: comment at the end of this document, Tweet at @xtimv, or email tim.f.vieira@gmail.com. If you found this article interesting please consider sharing it on social media. If you found this article useful please cite it

@misc{vieira2021order,
  author = {Tim Vieira},
  title = {On the Distribution Function of Order Statistics},
  year = {2021},
  url = {https://timvieira.github.io/blog/on-the-distribution-functions-of-order-statistics/}
}

Appendix¶

Symbolic tests for utility methods B and C

I have included some extra material for testing the efficient algorithms for summing over subsets of size $=K$ and size $\ge K$. This strategy compares the efficient and inefficient algorithms (brute-force enumeration) using symbolic mathematics. Using symbolic mathematics ensures that the tests do not depend on any specific numerical values. Furthermore, it is helpful during debugging as we can see which terms in the output expressions are different.

Inspecting the symbolic output

Below, we look at how the expression evaluated the subsets of size $=K$ in the B method. As a mnemonic, lower-case symbols represent the weight of inclusion into the set, and an upper-case symbols represent the weight of exclusion into the set. We consider a set with $3$ elements.

The different set sizes are the rows of the vector above. The first row is the weight of sets of size $0$, which is $X Y Z$ because all $3$ elements are excluded. The next line is we see that each term has one lower-case factor; this is because there is one element included in the set. See the pattern?

Finite-difference test

Below, we test gradients of the C function used in the pdf computation via the strategy described here.

Properties of $\nabla C$

Lemma: For any $j \in N$, the following differential identities hold: $$ \frac{\partial}{\partial p_j} C^N_K = C^{N \smallsetminus j}_{K-1} \quad\text{and}\quad \frac{\partial}{\partial q_j} C^N_K = C^{N \smallsetminus j}_{K} $$

We use these identities in our derivation of the PDF.