Sampling from the simplex

How do we sample uniformly at random from the unit dimensional simplex? It turns out there is an easy way to do it: generate i.i.d. uniform variables in [0,1] and let Y_{(k)} denote the k-th smallest among them. Then, (Y_{(1)}, Y_{(2)}-Y_{(1)}, Y_{(3)}-Y_{(2)}, \ldots, Y_{(n-1)}-Y_{(n-2)},1-Y_{(n-1)})  is a desired sample from the unit simplex \Delta_n=\{ (x_1,\ldots,x_n): x_i \geq 0, \sum_{i=1}^n x_i =1\}. It is straight-forward to prove now the claim and implement the method. However, there is a more efficient approach, based on the following exercise:

Exercise 1: Let Z_1,\ldots,Z_n be i.i.d. exponential random variables and define X_i = \frac{Y_i}{\sum_{j=1}^n Y_j}. Prove that (X_1,\ldots,X_n) is uniformly distributed in the (n-1) dimensional unit simplex \Delta_n.

As this paper explains in detail, the former method O(nlogn) time whereas the second O(n). Here is also a visualization of samples from the 2-simplex using a simple script.

Leave a comment