27
$\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$
7
  • 3
    $\begingroup$ This question may be relevant. $\endgroup$ Oct 8, 2013 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, 2013 at 11:42
  • 6
    $\begingroup$ Perhaps DirichletDistribution might help? $\endgroup$
    – chuy
    Oct 8, 2013 at 14:03
  • 1
    $\begingroup$ That question involved points on a sphere. Your constraint of $\sum{p_i}=1$ is different. $\endgroup$ Oct 8, 2013 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$ Feb 12, 2018 at 18:22

3 Answers 3

26
$\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$
6
  • $\begingroup$ Thanks a lot. Could you please explain in what respects does this behave differently from RandomVariate[UniformDistribution[]]? $\endgroup$
    – Schiphol
    Oct 9, 2013 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$ Oct 9, 2013 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, 2013 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$ Oct 11, 2013 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, 2013 at 18:28
11
$\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$
1
  • $\begingroup$ That's a slick capability. $\endgroup$
    – MikeY
    Jun 7, 2017 at 14:59
10
$\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$
1
  • 1
    $\begingroup$ 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! $\endgroup$ Aug 14, 2021 at 6:13

Your Answer

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

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