Re: Re: Re: Sample uniformly from a simplex
- To: mathgroup at smc.vnet.net
- Subject: [mg94554] Re: [mg94499] Re: [mg94473] Re: Sample uniformly from a simplex
- From: DrMajorBob <btreat1 at austin.rr.com>
- Date: Tue, 16 Dec 2008 02:34:19 -0500 (EST)
- References: <ghtjg8$rb3$1@smc.vnet.net> <200812141238.HAA10275@smc.vnet.net>
- Reply-to: drmajorbob at longhorns.com
Same 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
<jalex at lse.ac.uk> wrote:
> 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
>
--
DrMajorBob at longhorns.com
- References:
- Re: Sample uniformly from a simplex
- From: Jean-Marc Gulliet <jeanmarc.gulliet@gmail.com>
- Re: Sample uniformly from a simplex