Services & ResourcesWolfram ForumsMathGroup Archive

[Date Index] [Thread Index] [Author Index]

Sampling from a simplex: part 2

• To: mathgroup at smc.vnet.net
• Subject: [mg95816] Sampling from a simplex: part 2
• From: Andreas <aagas at ix.netcom.com>
• Date: Wed, 28 Jan 2009 06:28:39 -0500 (EST)

I posted a question back in Dec. about how to sample uniformly from a simplex.  (See: http://mathforum.org/kb/message.jspa?messageID=6530195&tstart=0 )

I received a bunch of interesting and helpful solutions.  Many thanks to all.

The solution I implemented looks like this (with a special hat tip to Dr. Major Bob):

sampleSimplex[sampleSize_, dimensions_] := Module[{s},
  s = RandomReal[ExponentialDistribution[2], {sampleSize, dimensions}];
  Table[Normalize[s[[i, All]], Total], {i, 1, sampleSize}]]

where: 
   sampleSize gives the number of points to sample and
   dimensions gives the number of coordinates to return for each point.

You can see the results graphically with this:

Graphics3D[Point /@ sampleSimplex[1000, 3], ImageSize -> 150, 
 BoxStyle -> Directive[LightGray], ViewPoint -> {Pi, Pi/2, 2}]

This works great for coordinate values between 0 and 1, but now I'd like to extend the approach beyond that range.

I came up with this:

sampleSimplexV2[sampleSize_, dimensions_, min_, max_] := 
 Module[{s, tableNormal},
  s = RandomReal[ExponentialDistribution[2], {sampleSize, dimensions}];
  tableNormal = 
   Table[Normalize[s[[i, All]], Total], {i, 1, sampleSize}];
  Table[min + tableNormal[[i]]*(max - min), {i, 1, sampleSize}]
  ] 

Where "min" and "max" set the range of the sides of the simplex.  

Example:

sampleSimplexV2[2, 3, -1, 2]

{{0.159696, 0.232216, -0.391911}, {-0.6432, 0.239791, 0.403409}}

Now my questions:

1. I think this works right.  Does it make sense?

2. Can I do this in a simpler, more direct, or more elegant way?

3. Question 2 opens up a larger question for me of approaches or rules of thumb to refactor Mathematica code (along the lines of Martin Fowler's book "Refactoring - Improving The Design of Existing Code" targeted predominately at object oriented languages) or simply to write the most concise, clear and robust code. Any suggestions, guidelines, or resource suggestions along these lines much appreciated.

A