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.

enter image description here

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. enter image description here


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$.

| cite | improve this answer | |
  • $\begingroup$ Gosh that's pretty slick -- and much better than anything I'd have come up with! $\endgroup$ – Daron Oct 13 '18 at 13:23
  • 1
    $\begingroup$ I've been trying to get some bounds for $X_k$ and hit the stage where I'm trying to show something like $5 a_1+ 4 a_2 + 3a_3 > a_4 + a_5$ with high probability, where $a_i$ are iid exponentials with mean $1$. As usual the PDF of the LHS is horrendous -- are there any good tricks to replace the LHS with something smaller and easier to work with? Maybe something like $12a$ for $a$ some exponential with mean $1$? $\endgroup$ – Daron Oct 13 '18 at 13:25
  • $\begingroup$ @Daron: I am not sure, because I am not an expert in this field. If I were you, I would ask a separate question about existence of tricks like that. $\endgroup$ – Kaban-5 Oct 13 '18 at 21:56
  • $\begingroup$ If you only need some bound (possibly quite loose) on the probability of stuff like $a_4 + a_5 > 5a_1 + 4a_2 + 3a_3$ happening, you may try calculating exponential moments. For example, let $t \in (0, 1)$ be some real number, then $\mathbb{P} (a_4 + a_5 > 5a_1 + 4a_2 + 3a_3) \le \mathbb{E} e^{t (a_4 + a_5 - 5a_1 - 4a_2 - 3a_3)} = \frac{1}{(1-t)^2 (1+3t)(1+4t)(1+5t)}$ and the last expression can be made less than $0.16$ by choosing $t$ carefully. But there are 2 problems with this approach: 1) bounds are quite crude 2) obtaining said bounds for all $n$ and $k$ at the time may be impossible. $\endgroup$ – Kaban-5 Oct 13 '18 at 22:06

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$.

| cite | improve this answer | |

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}$.

| cite | improve this answer | |

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service, privacy policy and cookie policy

Not the answer you're looking for? Browse other questions tagged or ask your own question.