MathGroup Archive 2008

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

Search the Archive

Re: Sample uniformly from a simplex

  • To: mathgroup at smc.vnet.net
  • Subject: [mg94473] Re: Sample uniformly from a simplex
  • From: Jean-Marc Gulliet <jeanmarc.gulliet at gmail.com>
  • Date: Sun, 14 Dec 2008 07:38:45 -0500 (EST)
  • Organization: The Open University, Milton Keynes, UK
  • References: <ghtjg8$rb3$1@smc.vnet.net>

Andreas wrote:

> 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

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



  • Prev by Date: Re: Sample uniformly from a simplex
  • Next by Date: Re: Import numerical data and split by pattern
  • Previous by thread: Re: Sample uniformly from a simplex
  • Next by thread: Re: Re: Sample uniformly from a simplex