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?

Improve this question
$\endgroup$
7
  • 3
    $\begingroup$ This question may be relevant. $\endgroup$ Commented 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$ Commented Oct 8, 2013 at 11:42
  • 6
    $\begingroup$ Perhaps DirichletDistribution might help? $\endgroup$ Commented 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$ Commented 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$ Commented Feb 12, 2018 at 18:22

3 Answers 3

Reset to default
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

Improve this answer
$\endgroup$
6
  • $\begingroup$ Thanks a lot. Could you please explain in what respects does this behave differently from RandomVariate[UniformDistribution[]]? $\endgroup$ Commented 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$ Commented 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$ Commented 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$ Commented 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$ Commented Oct 11, 2013 at 18:28
11
$\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

Improve this answer
$\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$ Commented Aug 14, 2021 at 6:13
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

Improve this answer
$\endgroup$
1
  • $\begingroup$ That's a slick capability. $\endgroup$ Commented Jun 7, 2017 at 14:59

Start asking to get answers

Find the answer to your question by asking.

Ask question

Explore related questions

See similar questions with these tags.