What is the probability distribution of the $k$th largest coordinate chosen over a simplex?

Question

Suppose we're selecting points uniformly at random from the $N$-simplex

$S_N = \{x \in \mathbb R^{N+1}: $ all $ x_i \ge 0$ and $x_1 + \ldots x_N = 1\}$.

One way to do this in practice is choose $N-1$ points $a_1, \ldots , a_{N-1}$ uniformly and independently from the unit interval $[0,1]$. Then for $a_0=0$ and $a_1=1$ construct the point $x \in S_N$ with each $x_i = a_i - a_{i-1}$.

One would expect for large $N$ the points $a_i$ to be evenly spaced across the interval, and so the the average point looks of $S_n$ looks pretty much like the constant $1/N$ vector. That is to say it's unlikely for any collection of entries to be small. Is anything known about the exact distribution of such collections?

Formally put suppose $X: \Omega \to S_N$ is a uniformly distributed random variable from some probability space onto the $N$-simplex. Define each $X_k : \Omega \to [0,1]$ by $X_k(x) = $ the $k$th largest coordinate of $X(x)$.

I've drawn some samples for $k = n/2$ which is the case I'm mostly interested in. It seems that, after you normalise the variable by multiplying by $N$, the mean tends to about $0.7$ (marked with a vertical line) from above. The distributions are also increasingly tighter bell-curves. The below is with 100,000 samples per curve.

Is there anything like a closed form known for the distribution of $X_k$? If not are there any useful bounds for probabilities like $P(X_k > 1/N \pm \epsilon)$ or $P(X_k < 1/N \pm \epsilon)$ that give answers similar to the behaviour above?

I am also interested in the variables $Y_k = X_1 + \ldots X_k$ if they are any easier to understand analytically, again primarily in the case $k = n/2$. Again the plots look like ever tighter bells and the mean tends to $0.15$ (the black line) from above.

Answer 1

Not sure if this is exactly what you need, but the following gives something that is close to the closed formula for distribution of $X_k$.

As noted by Mark Meckes, the random point $x$ from the simplex has the same distribution as $\left( \dfrac{b_1}{b_1 + \ldots + b_n}, \dfrac{b_2}{b_1 + \ldots + b_n}, \ldots, \dfrac{b_n}{b_1 + \ldots + b_n} \right )$, where $b_i$ are independent exponential random variables with expectation $1$.

Now, we want to know the distribution of $X_k$, the $k$-th largest coordinate of $x$. Another well known result is that vector $(B_n, B_{n - 1}, \ldots B_1)$, where $B_k$ is the $k$-th largest coordinate of $(b_1, b_2, \ldots b_n)$ (recall that $b_i$ are independent exponential variables with expectation $1$) has the same distribution as $ \left( \dfrac{\xi_n}{n}, \dfrac{\xi_n}{n} + \dfrac{\xi_{n - 1}}{n - 1}, \ldots, \dfrac{\xi_n}{n} + \dfrac{\xi_{n-1}}{n-1} + \ldots + \dfrac{\xi_1}{1} \right )$, where $\xi_1, \xi_2, \ldots, \xi_n$ are independent exponential random variables with expectation $1$. (You can find the proof here, for example).

It follows that $X_k$ has the distribution $\dfrac{\xi_n / n + \xi_{n-1} / (n-1) + \ldots + \xi_k / k}{\xi_n + \xi_{n - 1} + \ldots + \xi_1}$. This expression is still not pretty, but it allows us to find, for example, $\mathbb{E} X_k$. To find $\mathbb{E} X_k$, notice that independence of $\xi_i$ implies that all random variables of form $\dfrac{\xi_{p_n} / n + \xi_{p_{n-1}} / (n-1) + \ldots + \xi_{p_k} / k}{\xi_n + \xi_{n - 1} + \ldots + \xi_1}$, where $p$ is a permutation of integers from $1$ to $n$, have the same expectation. Therefore, $$\mathbb{E} X_k = \mathbb{E} \dfrac{\xi_n / n + \xi_{n-1} / (n-1) + \ldots + \xi_k / k}{\xi_n + \xi_{n - 1} + \ldots + \xi_1} = \dfrac{1}{n!} \times\sum\limits_{p \in S_n} \mathbb{E} \dfrac{\xi_{p_n} / n + \xi_{p_{n-1}} / (n-1) + \ldots + \xi_{p_k} / k}{\xi_n + \xi_{n - 1} + \ldots + \xi_1} = \dfrac{1}{n} \left( \dfrac{1}{n} + \ldots + \dfrac{1}{k} \right) \dfrac{\xi_1 + \xi_2 + \ldots + \xi_n}{\xi_1 + \xi_2 + \ldots + \xi_n} = \frac{1}{n} (H_n - H_{k - 1}),$$ where $H_m := 1 + \frac{1}{2} + \ldots + \frac{1}{m}$. Here we used averaging over all permutations to make coefficients before each $\xi_i$ equal to each other.

So, $\lim\limits_{n \to +\infty} \mathbb{E} \frac{X_{n}}{2n} = \lim\limits_{n \to +\infty} (H_{2n} - H_{n - 1}) = \lim \limits_{n \to +\infty} (\ln (2n) - \ln (n - 1)) = \ln 2 \approx 0.7$, which explains the first phenomena you observed. Writing down $\mathbb{E} \frac{Y_n}{2n}$ in the same way should explain the second phenomena you observed.

Now, I am not sure how good is the representation above for your purposes, but studying them already can give some interesting results about distributions of $X_k$.

Answer 2

A little too long for a comment, but I don't have time right now to turn this observation into a proper answer:

Another way to generate a uniform random point in the simplex is to let $b_1, \ldots, b_n$ be independent exponential random variables, and let $$ x_i = \frac{b_i}{b_1 + \cdots + b_n}. $$ It then turns out that the point $(x_1, \ldots, x_n)$ in the simplex is independent of the Gamma-distributed random variable $z_n = b_1 + \cdots + b_n$, and of course if you define $B_k$ to be the $k$th largest $b_i$, then $X_k = B_k/Z_n$.

Answer 3

As proposed before, we can write $x_i$ as $$x_i=\frac{b_i}{b_1+\cdots+b_n} $$ where $b_i$ are independent exponential variable with parameter one.

Let $y_i=e^{-b_i}$ then $y_i$ are iid uniform random variable on $[0,1]$. Then the law of the $k$th larger $y_i$ follow beta law $B(n-k+1,k)$ https://en.wikipedia.org/wiki/Beta_distribution#Order_statistics. (Indeed, because there are $k-1$ elements in $[y,1]$ and $n-k$ elements in $[0,y]$ the density is proportionnal to $y^{n-k-1}(1-y)^{k}$ )

For large $n$ as $\frac{1}{n}\sum b_i\rightarrow 1$ in probability then the law of $n x_i$ converge to $-\log (y_i)$. And for example $X_{\frac{n}{2}}$ has mean $\log(2)$ has $Y_{\frac{n}{2}}\rightarrow \frac{1}{2}$.

And then you can use all the results from beta distribution. (I think it is the best for your "anything like a closed form")

For your last question, once you fixe $Y_k$, all the $(Y_i)_{i> k}$ are uniform on $[0,Y_k]$. And then the mean should be $\sum_{i=1}^k X_i=-\sum_{i=1}^k \log(Y_i) \sim -k $ $\int_0^{Y_k} \log(x)dx$. The integral gives for $\frac{1-\log(2)}{2}\approx 0.153 $ for $k=\frac{n}{2}$.

Answer 4

It turns out that Kaban-5's answer (which is excellent) can be extended into analytical density for the lowest element of a uniformly distributed simplex, which we can extend to densities of $j$-th lowest element given all the smaller elements. This turns out to be enough for some statistical applications of ordered simplices, so sharing here if anybody besides me finds it useful.

We'll start with Kaban-5's result:

$$ X_k = \dfrac{\xi_n / n + \xi_{n-1} / (n-1) + \ldots + \xi_k / k}{\xi_n + \xi_{n - 1} + \ldots + \xi_1} $$

where $\xi_i \sim \mathrm{Exponential}(1)$

For $X_n$, the lowest element this simplifies to:

$$ X_n = \dfrac{\xi_n}{n(\xi_n + \xi_{n - 1} + \ldots + \xi_1)} = \dfrac{\xi_n}{n(\xi_n + \gamma)} $$

where $\xi_n$ and $\gamma$ are independent and $\gamma \sim \mathrm{Gamma(n - 1, 1)}$.

We then have

$$ Pr(X_n < x) = Pr(\dfrac{\xi_n}{n(\xi_n + \gamma)} < x) = Pr(\gamma > \frac{\xi_n - \xi_n n x}{nx}) = \\ = \int_0^\infty f_{\xi_n}(x) \left(1 - F_\gamma \left( \frac{\xi_n - \xi_n n x}{nx}\right)\right) \mathrm{d}x $$

where $f_{\xi_n}$ is the PDF of $\xi_n$ (exponential distribution) and $F_\gamma$ is the CDF of $\gamma$. My good friend Wolfram Alpha tells me that under such circumstances we get CDF for $X_n$ as

$$ Pr(X_n < x) = \frac{nx - 1 + (1 - nx)^n}{nx - 1} $$

which we can differenatiate to get the PDF as $f_{X_n}(x) = n(n - 1) (1 - nx)^{n - 2} : 0 < x < \frac{1}{n}$.

We further observe that if $X_n, ..., X_1$ is a uniformly distributed ordered simplex with smallest element $X_n$, then $\frac{X_{n-1} - X_n}{1 - nX_n}, ..., \frac{X_{1} - X_n}{1 - nX_n}$ is also a uniformly distributed ordered simplex of size $n - 1$

This is actually sufficient for (some) applications: we can now recursively build the marginal density of $j$-th smallest element given all the smaller ones, which is just enough to use ordered simplex as a building block in statistical models fitted with (Hamiltonian) MCMC (which is what I'm interested in). See https://discourse.mc-stan.org/t/ordered-simplex-constraint-transform/24102/14 for an implementation in the Stan probabilistic programming language.