mathematica, Sample uniformly from a simplex

I need to develop Mathematica code to sample uniformly from a unit n-dimensional simplex.

I came across a description of the problem at: http://geomblog.blogspot.com/2005/10/sampling-from-simplex.html

Specifically, I would like a uniform sample from the set

X = { (x1, x2, ..., xD) | 0 <= xi <= 1, x1 + x2 + ... + xD = 1}.

D is the dimension of the simplex.

So, the coordinates of any point on the simplex would sum to 1 and I need to sample points on the simplex.

geomblog's solution suggested:

Generating IID random samples from an exponential distribution by sampling X from [0,1] uniformly, and returning -log(X)).

Take n samples, then normalize.

This should result in a list of numbers which is a uniform sample from the simplex.

I've searched extensively for a Mathematica implementation of something like this, to no avail.

I keep trying different things but haven't made much headway.

Any suggestions for how to develop this (or an equivelant) in Mathematica much appreciated

A

I've kept thinking about this and came up with the following.

sampleSimplex[sampleSize_, dimensions_] :=
Module[{sample, sampleOutput},
sample = RandomReal[1, {sampleSize, dimensions}];
sampleOutput =
Table[1 / Total[sample[[i, All]]] *sample[[i, All]], {i, 1,
sampleSize}];
sampleOutput
]

example:

sampleSimplex[4, 3]

{{0.163972, 0.445817, 0.390211}, {0.470828, 0.291766,
0.237406}, {0.56663, 0.3813, 0.0520696}, {0.182715, 0.600876, 0.216409}}

This doesn't draw from an exponential distribution, but would this work?

Also, if it does work can I do it a better (faster, more elegant) way?

nsional simplex.
m/2005/10/sampling-from-simplex.html
to sample points on the simplex.
g X from [0,1] uniformly, and returning -log(X)).
e simplex.
ike this, to no avail.
much appreciated

Hi. Try this (using Version 6 or 7):

ranSimp[d_] := (#/Plus @@ #) &[-Log[RandomReal[1, d]]]
ranSimp[d_, n_] := (#/Plus @@ #) & /@ (-Log[RandomReal[1, {n, d}]])

d is the "dimension" of the simplex and n is the number of draws.

--Mark

nsional simplex.
m/2005/10/sampling-from-simplex.html
to sample points on the simplex.
g X from [0,1] uniformly, and returning -log(X)).
e simplex.
ike this, to no avail.
much appreciated

The following is code from Mathematica tutorials for sampling from the
Dirichlet distribution.

DirichletDistribution /: Random`DistributionVector[
DirichletDistribution[alpha_?(VectorQ[#, Positive] &)], n_Integer,
prec_?Positive] :=
Block[{gammas},
gammas =
Map[RandomReal[GammaDistribution[#, 1], n,
WorkingPrecision -> prec] &, alpha];
Transpose[gammas]/Total[gammas]]

Asim

The following example in R^3 should be close to what you are looking
for, though I am not a statistician and I may have failed to fully grasp
the suggested solution.

In[1]:= Module[{x},
Table[x = -Log[RandomReal[1, {3}]]; x/Total[x], {5}]]
Total[Transpose[%]]

Out[1]= {{0.545974, 0.204439, 0.249587}, {0.36947, 0.0545329,
0.575997}, {0.523704, 0.319784, 0.156512}, {0.490651, 0.398176,
0.111173}, {0.0332044, 0.406806, 0.55999}}

Out[2]= {1., 1., 1., 1., 1.}

In[3]:= Graphics3D[
Point /@ Module[{x},
Table[x = -Log[RandomReal[1, {3}]]; x/Total[x], {1000}]]]

Hope this helps,
-- Jean-Marc

Yet another solution would be to use a "stick-breaking" algorithm.
Generate N-1 random reals in the interval [0,1], then use the length
of the N subintervals as the values.

randomSimplex[n_] := (#[[2]] - #[[1]]) & /@ (Partition[Union[{0, 1},
RandomReal[{0, 1}, n - 1]], 2, 1])

You can also use this method to sample points from a convex polytope
with a little extra work. See http://cg.scs.carleton.ca/ ~luc/rnbookindex.html
, page 568 (theorem 2.1).

Cheers,

Jason










Please access the attached hyperlink for an important electronic communications disclaimer: http://www.lse.ac.uk/collections/secretariat/legal/disclaimer.htm

ame idea slightly simpler:

randomSimplex[n_] :=
Rest@# - Most@# &[Union[{0, 1}, RandomReal[{0, 1}, n - 1]]]

randomSimplex[10]

{0.00288906, 0.00240762, 0.0739569, 0.252289, 0.188723, 0.00997999, \
0.10469, 0.172935, 0.114573, 0.0775569}

Bobby

On Mon, 15 Dec 2008 06:43:14 -0600, J. McKenzie Alexander
< XXXX@XXXXX.COM > wrote:




--
XXXX@XXXXX.COM

Great set of solutions. Many thanks to all who contributed.
A