Uniformly distributed n-dimensional probability vectors over a simplex

Question

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?

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 08 '13 at 11:42
That question involved points on a sphere. Your constraint of $\sum{p_i}=1$ is different. — Sjoerd C. de Vries, Oct 08 '13 at 14:39
I agree with chuy. From wikipedia: Dirichlet Distribution, we have: "When alpha = 1, the symmetric Dirichlet distribution is equivalent to a uniform distribution over the open standard K-1-simplex, i.e. it is uniform over all points in its support". Using DirichletDistribution is most likely then also the best implementation for most purposes. — Jacob Akkerboom, Oct 08 '13 at 15:26
Related: Generate regularly spaced points from the surface of simplexes. — István Zachar, Oct 11 '13 at 14:16
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. https://stackoverflow.com/questions/18659858/generating-a-list-of-random-numbers-summing-to-1 — cgnorthcutt, Feb 12 '18 at 18:22

Ray Koopman · Accepted Answer · 2013-10-14T21:02:19.417

26

#/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[[2]]-p[[1]], Sqrt@3*p[[3]]} are the barycentric coordinates of the 3-element probability vector p:

enter image description here

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:

enter image description here

edited Oct 14 '13 at 21:02

answered Oct 08 '13 at 15:49

Ray Koopman

3,306
14
13

Thanks a lot. Could you please explain in what respects does this behave differently from RandomVariate[UniformDistribution[]]? – Schiphol Oct 09 '13 at 08:10
See for yourself. Try it with n = 2 and make a histogram of p[[1]]. Or use n = 3 and ListPlot the barycentric coordinates: {x,y} = {p[[2]]-p[[1]],Sqrt@3*p[[3]]}. – Ray Koopman Oct 09 '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
5

You can also use Mathematica's built-in DirichletDistribution: points = RandomVariate[DirichletDistribution[{1, 1, 1}], 2000] /. v_?VectorQ :> {v[[2]] - v[[1]], Sqrt[3] (1 - Total[v])}; and then ListPlot[points]. – chuy Oct 11 '13 at 18:28
Amazing! Does anyone have a prove that this distribution is uniform on the simplex? It passes every stat sampling test I tried, but I still can't prove it's right. I posted this question on the Math page too. https://math.stackexchange.com/questions/4235192/prove-this-method-generates-uniform-points-on-the-simplex – Jerry Guern Aug 28 '21 at 17:22

score 11 · Answer 2 · answered Jun 06 '17 at 21:26

11

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]]

answered Jun 06 '17 at 21:26

Carl Woll

130,679
6
243
355

That's a slick capability. – MikeY Jun 07 '17 at 14:59

MikeY · Answer 3 · 2017-06-06T19:50:50.450

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…

Generate d-1 uniformly distributed random values in the range [0,1]
Add a 0 and a 1 to the list
Sort the list
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

This is more or less what I came up with on my own, but then I found this Question page while looking for a more efficient way! — Jerry Guern, Aug 14 '21 at 06:13

Uniformly distributed n-dimensional probability vectors over a simplex

3 Answers3

Linked

Related