# Uniform sampling from a simplex

I am looking for an algorithm to generate an array of N random numbers, such that the sum of the N numbers is 1, and all numbers lie within 0 and 1. For example, N=3, the random point (x, y, z) should lie within the triangle:

x + y + z = 1
0 < x < 1
0 < y < 1
0 < z < 1


Ideally I want each point within the area to have equal probability. If it's too hard, I can drop the requirement. Thanks.

• What is the target distribution? What have you tried? – Raphael Aug 16 '12 at 22:07
• Note that there is always rejection sampling: sample $n$ uniform numbers and reject if the numbers don't add up to $1$. Here, the expected number of iterations is uncomfortably high, so you should do something else. – Raphael Aug 16 '12 at 23:14

Let us first assume that you want to sample within

x + y + z = 1
0 ≤ x ≤ 1
0 ≤ y ≤ 1
0 ≤ z ≤ 1


This doesn't make quite a difference, since the sample point will still lie in your requested area with high probability.

Now you are left with sampling a point from a simplex. In the 3d example you get a 2d simplex (triangle) realized in 3d.

How to pick a point uniformly at random was discussed in this blog post (see the comments).

For your problem it would mean that you take $n-1$ random numbers from the interval $(0,1)$, then you add a $0$ and $1$ to get a list of $n+1$ numbers. You sort the list and then you record the differences between two consecutive elements. This gives you a list of $n$ number that will sum up to $1$. Moreover this sampling is uniform. This idea can be found in Donald B. Rubin, The Bayesian bootstrap Ann. Statist. 9, 1981, 130-134.

For example ($n=4$) you have the three random numbers 0.4 0.2 0.1 then you obtain the sorted sequence 0 0.1 0.2 0.4 1 and this gives the differences 0.1 0.1 0.2 0.6, and by construction these four numbers sum up to 1.

Another approach is the following: first sample from the hypercube (that is you forget about x+y+z=1) and then normalize the sample point. The normalization is a projection from the $d$-hypercube to the $d-1$-simplex. It should be intuitively clear that the points at the center of the simplex have more "pre-image-points" than at the outside. Hence, if you sample uniformly from the hypercube, this wont give you a uniform sampling in the simplex. However, if you sample from the hypercube with an appropriate Exponential Distribution, than this effect cancels out. The Figure gives you an idea how both methods will sample. However, I prefer the "sorting" method due to its simple form. It's also easier to implement.

• I guess the naive idea -- draw $n$ numbers from $(0,1)$ and normalise -- is faulty, then. – Raphael Aug 16 '12 at 23:13
• Are there a simple proof that shows the sorting one gives a uniform distribution? I have only elementary background in probability so the paper is over my head. – Chao Xu Aug 17 '12 at 21:52
• @ChaoXu just observe that there is a 1-to-1 correspondence between $n$ numbers in the simplex and partitions of the interval $(0, 1)$ into $n$ subintervals. the sampling algorithm corresponds to "throwing" $n-1$ random "interval endpoints" on $(0, 1)$. you can try to verify that that's uniform, say by induction and using conditional probability. – Sasho Nikolov Aug 19 '12 at 19:21
• @Orient: Please ask you questions in a separate post and dont misuse the comments for this. – A.Schulz Jul 17 '14 at 11:36

Devroye is an excellent reference for questions of this sort. Chap.7 gives the algorithms needed to generate uniform order statistics, which the OP is after.

For generating uniform order statistics, sorting $n$ samples of $[0,1]$ will do. This approach takes $O(n \log n)$ time. A quicker way ( available in the book) involves sampling $n$ random numbers $x_1,\ldots,x_n$ from an $\mathrm{Exp}(1)$ pdf. (These are the spacings of the uniform pdf). Then, return the values $$(y_i)_{1\leq i\leq n} = \frac{\sum \limits_{1\ldots i} x_j}{\sum \limits_{1\ldots n} x_j}$$ which are automatically sorted, in $O(n)$ time overall. (I'm overlapping with A.Schulz's answer here- just making the computation more explicit).

The same approach can be adapted, via inverse CDF Sampling, to sample any non-uniform pdf over $[0,1]$. There's also a trick that enables you to sample uniformly over a simplex other than the canonical simplex (say $2x+3y+z = 5$).

• If I follow the answer here: stackoverflow.com/questions/2106503/… Then generating random number from exponential distribution involves evaluating logarithm, which can be a bit slow. – R zu Nov 7 '18 at 19:32
X[0] = 0
for i = 1 to N-1
X[i] = uniform(0,1)
X[n] = 1
sort X[0..N]
for i = 1 to N
Z[i] = X[i] - X[i-1]
return Z[1..N]


Here, uniform(0,1) returns a real number independently and uniformly distributed between 0 and 1.

• This is A. Schulz's answer in code without the explanation, right? – Raphael Sep 28 '12 at 14:21

See this paper: Smith, N. and Tromble, R., Sampling uniformly from the unit simplex.

• Please format your answer in a readable way: you're writing for human beings, not the bibtex compiler. Also, if the paper's available online, it's much more efficient for you to provide a link. – David Richerby Dec 9 '15 at 7:43