24
$\begingroup$

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?

$\endgroup$
  • 3
    $\begingroup$ This question may be relevant. $\endgroup$ – Sjoerd C. de Vries Oct 8 '13 at 11:21
  • $\begingroup$ 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. $\endgroup$ – Schiphol Oct 8 '13 at 11:42
  • 6
    $\begingroup$ Perhaps DirichletDistribution might help? $\endgroup$ – chuy Oct 8 '13 at 14:03
  • 1
    $\begingroup$ That question involved points on a sphere. Your constraint of $\sum{p_i}=1$ is different. $\endgroup$ – Sjoerd C. de Vries Oct 8 '13 at 14:39
  • 2
    $\begingroup$ 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/… $\endgroup$ – cgnorthcutt Feb 12 '18 at 18:22
24
$\begingroup$

#/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

$\endgroup$
  • $\begingroup$ Thanks a lot. Could you please explain in what respects does this behave differently from RandomVariate[UniformDistribution[]]? $\endgroup$ – Schiphol Oct 9 '13 at 8:10
  • $\begingroup$ 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]]}. $\endgroup$ – Ray Koopman Oct 9 '13 at 18:41
  • $\begingroup$ 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. $\endgroup$ – Schiphol Oct 11 '13 at 10:39
  • $\begingroup$ 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. $\endgroup$ – Ray Koopman Oct 11 '13 at 13:55
  • 5
    $\begingroup$ 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]. $\endgroup$ – chuy Oct 11 '13 at 18:28
7
$\begingroup$

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.

enter image description here

This topic is well-covered here... https://stackoverflow.com/questions/3010837/sample-uniformly-at-random-from-an-n-dimensional-unit-simplex

$\endgroup$
6
$\begingroup$

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

enter image description here

$\endgroup$
  • $\begingroup$ That's a slick capability. $\endgroup$ – MikeY Jun 7 '17 at 14:59

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service, privacy policy and cookie policy

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