# Uniformly distributed n-dimensional probability vectors over a simplex

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?

-
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
Perhaps DirichletDistribution might help? –  chuy Oct 8 '13 at 14:03
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

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