Re: Re: Sample uniformly from a simplex
- To: mathgroup at smc.vnet.net
- Subject: [mg94499] Re: [mg94473] Re: Sample uniformly from a simplex
- From: "J. McKenzie Alexander" <jalex at lse.ac.uk>
- Date: Mon, 15 Dec 2008 07:43:14 -0500 (EST)
- References: <ghtjg8$rb3$1@smc.vnet.net> <200812141238.HAA10275@smc.vnet.net>
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
--
Dr J. McKenzie Alexander
Department of Philosophy, Logic and Scientific Method
London School of Economics and Political Science
Houghton Street, London WC2A 2AE
On 14 Dec 2008, at 12:38, Jean-Marc Gulliet wrote:
> 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
>
>
Please access the attached hyperlink for an important electronic communications disclaimer: http://www.lse.ac.uk/collections/secretariat/legal/disclaimer.htm
- References:
- Re: Sample uniformly from a simplex
- From: Jean-Marc Gulliet <jeanmarc.gulliet@gmail.com>
- Re: Sample uniformly from a simplex