Take the 2-minute tour ×
Mathematica Stack Exchange is a question and answer site for users of Mathematica. It's 100% free, no registration required.

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?

share|improve this question
3  
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
3  
Perhaps DirichletDistribution might help? –  chuy Oct 8 '13 at 14:03
1  
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
    
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 8 '13 at 15:26
show 1 more comment

1 Answer

up vote 8 down vote accepted

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

share|improve this answer
    
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[[1]]. Or use n = 3 and ListPlot the barycentric coordinates: {x,y} = {p[[2]]-p[[1]],Sqrt@3*p[[3]]}. –  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
1  
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
add comment

Your Answer

 
discard

By posting your answer, you agree to the privacy policy and terms of service.

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