36
$\begingroup$

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.

Improve this question
$\endgroup$
2
  • $\begingroup$ What is the target distribution? What have you tried? $\endgroup$
    – Raphael
    Commented Aug 16, 2012 at 22:07
  • 3
    $\begingroup$ 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. $\endgroup$
    – Raphael
    Commented Aug 16, 2012 at 23:14

6 Answers 6

Reset to default
37
$\begingroup$

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.

Example of the 2 sampling methods

Improve this answer
$\endgroup$
17
  • $\begingroup$ I guess the naive idea -- draw $n$ numbers from $(0,1)$ and normalise -- is faulty, then. $\endgroup$
    – Raphael
    Commented Aug 16, 2012 at 23:13
  • $\begingroup$ I addressed your question in the extended answer. $\endgroup$
    – A.Schulz
    Commented Aug 17, 2012 at 8:59
  • 1
    $\begingroup$ 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. $\endgroup$
    – Chao Xu
    Commented Aug 17, 2012 at 21:52
  • 5
    $\begingroup$ @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. $\endgroup$
    – Sasho Nikolov
    Commented Aug 19, 2012 at 19:21
  • 1
    $\begingroup$ @Orient: Please ask you questions in a separate post and dont misuse the comments for this. $\endgroup$
    – A.Schulz
    Commented Jul 17, 2014 at 11:36
9
$\begingroup$

This is to add to the existing answers.

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

Improve this answer
$\endgroup$
2
  • $\begingroup$ 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. $\endgroup$
    – R zu
    Commented Nov 7, 2018 at 19:32
  • $\begingroup$ Note sure if its a mistake, but I believe the limits should read $(y_i)_{1\leq i < n}$, so that we get $n-1$ sorted Uniform samples [0,1] from $n$ Exponential(1) samples? If not, it seems to me, that the last sample would always equal 1. $\endgroup$
    – chrish.
    Commented Oct 28, 2022 at 18:14
3
$\begingroup$ 
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.

Improve this answer
$\endgroup$
1
  • 6
    $\begingroup$ This is A. Schulz's answer in code without the explanation, right? $\endgroup$
    – Raphael
    Commented Sep 28, 2012 at 14:21
2
$\begingroup$

One other possibility is to use the Dirichlet distribution using the SciPy module. Essentially, if we set $\alpha=(1,\dots,1)^\top$, then the probability density is effectively uniform on a $n$-dimensional simplex where $n$ is the number of elements in $\alpha$.

Here is an example:

import numpy as np
from scipy.stats import dirichlet
n = 2
size = 1000
alpha = np.ones(n)
samples = dirichlet.rvs(size=size, alpha=alpha)

Now, if I plot it using matplotlib, we can see immediately see that it does a good job:

import matplotlib.pyplot as plt
plt.ion()
plt.scatter(samples[:,0], samples[:,1], alpha=0.05)

An illustration of sampling from a two-dimensional Dirichlet distribution.

Improve this answer
$\endgroup$
1
$\begingroup$

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

Improve this answer
$\endgroup$
1
  • 2
    $\begingroup$ 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. $\endgroup$
    – David Richerby
    Commented Dec 9, 2015 at 7:43
0
$\begingroup$ 
import numpy as np
def rand_simplex(k):
    return tuple(np.random.dirichlet((1,)*k))

This will generate one sample at a time.

Improve this answer
$\endgroup$
2
  • 1
    $\begingroup$ Could you please explain the algorithm behind np.random.dirichlet? $\endgroup$
    – xskxzr
    Commented Dec 10, 2020 at 2:50
  • 1
    $\begingroup$ It produces samples from a dirichlet distribution en.wikipedia.org/wiki/Dirichlet_distribution. Intuitively, it is "stick breaking": Take a stick of unit length and randomly break it resulting in k pieces. The lengths of those pieces is the random vector. The parameters are "concentrations" governing the tendencies of each element of the vector to vary between samples. $\endgroup$
    – Harris
    Commented Dec 14, 2020 at 18:05

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service and acknowledge you have read our privacy policy.

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