Uniformly distributed n-dimensional probability vectors over a simplex

What's the right way to generate a random probability vector $p={p_1,\ldots,p_n} \in {(0,1)}^n$ where $\sum_i p_i=1$, uniformly distributed over the $(n-1)$-dimensional simplex?

What I have is

Intervals = Table[{0, 1}, {i, n}]
RandomPoint := Block[{a},
a = RandomVariate[UniformDistribution[Intervals]];
a/Total[a]];

But I am unsure that this is correct. In particular, I'm unsure that it's any different from:

RandomPoint := Block[{a},
a = Table[Random[], {i, n}];
a/Total[a]];

And the latter clearly will not distribute vectors uniformly. Is the first code the right one?

• This question may be relevant. – Sjoerd C. de Vries Oct 8 '13 at 11:21
• Thanks, @SjoerdC.deVries. That question seems to suggest that my first code is also incorrect? I'm assuming that that bunch of smart guys would have stumbled upon it. – Schiphol Oct 8 '13 at 11:42
• Perhaps DirichletDistribution might help? – chuy Oct 8 '13 at 14:03
• That question involved points on a sphere. Your constraint of $\sum{p_i}=1$ is different. – Sjoerd C. de Vries Oct 8 '13 at 14:39
• Some folks here might find this StackOverflow duplicate useful (which is asking exactly the same question, but from a computer science perspective). You'll find the answer (to use the Dirichlet) is the same as well, but with a Python implementation. stackoverflow.com/questions/18659858/… – cgnorthcutt Feb 12 '18 at 18:22

#/Total[#,{2}]&@Log@RandomReal[{0,1},{m,n}] will give you a sample of m points from a uniform distribution over an n-1-dimensional regular simplex. (An equilateral triangle is a 2-dimensional regular simplex.) Here's what m = 2000, n = 3 should look like, where {x,y} = {p[]-p[], Sqrt@3*p[]} are the barycentric coordinates of the 3-element probability vector p: Here's what you get if you omit the Log@ and normalize Uniform(0,1) variables, which is what both of the OP's examples do: • Thanks a lot. Could you please explain in what respects does this behave differently from RandomVariate[UniformDistribution[]]? – Schiphol Oct 9 '13 at 8:10
• See for yourself. Try it with n = 2 and make a histogram of p[]. Or use n = 3 and ListPlot the barycentric coordinates: {x,y} = {p[]-p[],Sqrt@3*p[]}. – Ray Koopman Oct 9 '13 at 18:41
• Yes, the difference is clear -- see my answer below. Actually, I meant for you to explain the difference in algorithmic terms, or perhaps provide pointers to a textbook explanation of why your method is doing what it's doing. – Schiphol Oct 11 '13 at 10:39
• I generate a Dirichlet distribution in which all the concentration parameters are 1. See the link that Jacob provided, then scroll down to this section and remember that the log of a Uniform(0,1) variable is proportional to a Gamma variable with shape parameter 1. – Ray Koopman Oct 11 '13 at 13:55
• You can also use Mathematica's built-in DirichletDistribution: points = RandomVariate[DirichletDistribution[{1, 1, 1}], 2000] /. v_?VectorQ :> {v[] - v[], Sqrt (1 - Total[v])}; and then ListPlot[points]. – chuy Oct 11 '13 at 18:28

Old question, but I didn't see this method. Generates $n$ points uniformly randomly distributed on a simplex embedded in $d$ dimensions.

genSimplex[n_, d_] :=
Table[Differences[Sort[Flatten[{0, RandomReal[1, d – 1], 1}]]], {n}];

The algorithm generates points that are randomly distributed on an outer face of a simplex. The way to generate them is, for a d-dimensional problem…

1. Generate d-1 uniformly distributed random values in the range [0,1]
2. Add a 0 and a 1 to the list
3. Sort the list
4. Extract a list of the differences between the elements

You now have a list of random values that sum to 1 (so they are on a plane that is defined by points that sum to one) and that are otherwise independent of each other, so their dispersion is uniform.

Updating the answer with a picture of example data with 1,000 points. This topic is well-covered here... https://stackoverflow.com/questions/3010837/sample-uniformly-at-random-from-an-n-dimensional-unit-simplex

Starting in M10.2, you can just use RandomPoint:

pts=RandomPoint[Simplex[{{0,0,1},{0,1,0},{1,0,0}}], 1000];
Graphics3D[Point[pts]] • That's a slick capability. – MikeY Jun 7 '17 at 14:59