Path: news.mathworks.com!not-for-mail
From: <HIDDEN>
Newsgroups: comp.soft-sys.matlab
Subject: Re: Generate random numbers with fixed sum and different constraints
Date: Tue, 20 Nov 2012 19:54:09 +0000 (UTC)
Organization: The MathWorks, Inc.
Lines: 43
Message-ID: <k8gn51$42$1@newscl01ah.mathworks.com>
References: <k80995$9eo$1@newscl01ah.mathworks.com> <k87s11$h6u$1@newscl01ah.mathworks.com> <k88l7g$7tm$1@newscl01ah.mathworks.com> <k88qag$nm1$1@newscl01ah.mathworks.com> <k89ehg$rm6$1@newscl01ah.mathworks.com> <k8b4e9$ja7$1@newscl01ah.mathworks.com> <k8badi$8rd$1@newscl01ah.mathworks.com> <k8ck6h$lhp$1@newscl01ah.mathworks.com> <k8d3bv$a2s$1@newscl01ah.mathworks.com> <k8d77s$lfd$1@newscl01ah.mathworks.com> <k8d8mk$pv7$1@newscl01ah.mathworks.com> <k8d99h$rrh$1@newscl01ah.mathworks.com> <k8dakj$39u$1@newscl01ah.mathworks.com> <k8ee8q$k6f$1@newscl01ah.mathworks.com> <k8f9ck$gng$1@newscl01ah.mathworks.com>
Reply-To: <HIDDEN>
NNTP-Posting-Host: www-03-blr.mathworks.com
Content-Type: text/plain; charset=UTF-8; format=flowed
Content-Transfer-Encoding: 8bit
X-Trace: newscl01ah.mathworks.com 1353441249 130 172.30.248.48 (20 Nov 2012 19:54:09 GMT)
X-Complaints-To: news@mathworks.com
NNTP-Posting-Date: Tue, 20 Nov 2012 19:54:09 +0000 (UTC)
X-Newsreader: MATLAB Central Newsreader 1187260
Xref: news.mathworks.com comp.soft-sys.matlab:783170

"Bruno Luong" <b.luong@fogale.findmycountry> wrote in message <k8f9ck$gng$1@newscl01ah.mathworks.com>...
> Really cool Roger! I would love if you could share the justification of this formula.

"james bejon" wrote in message <k8fe8e$2lo$1@newscl01ah.mathworks.com>...
> Very neat!  (Though I have no idea how it works).
- - - - - - - -
  To Bruno, James, and anyone else interested:

  I recommended the two lines of code

 R = bsxfun(@power,rand(m-1,n),1./(m-1:-1:1)');
 X = cumprod([ones(1,n);R]).*[ones(m,n)-[R;zeros(1,n)]];

as a substitute for Bruno's

 w = randfixedsum(m,n,1,0,1)

(with X here instead of w) because 'randfixedsum' is designed for handling complicated polytopes which break down to into a great many individual simplices, whereas Bruno's call represents a single m-1 dimensional simplex within R^m.  (It's analogous to using a cannon to swat flies.)  It is more efficient to compute the single simplex directly without having to bother with provisions for a decomposition into simplices.  Its vertices are the m points:

 [1;0;0;0;...], [0;1;0;0;...], [0;0;1;0;...], ..., [0;0;0;...;0;0;1]

which we can call P1, P2, ..., Pm.

  For example letting m = 4 we have a 3D tetrahedron simplex imbedded in R^4 space, and as can be seen the two code lines, in effect, carry out the computation

 r1 = rand^(1/3); r2 = rand^(1/2); r3 = rand^(1);
 X = [ (1)*(1-r1) ; (1*r1)*(1-r2) ; (1*r1*r2)*(1-r3) ; (1*r1*r2*r3)*(1-0) ];

which with the definitions of the points is the same as

 X = (1-r1)*P1 + r1*((1-r2)*P2 + r2*((1-r3)*P3 + r3*P4));

  This distributes points within the tetrahedron in a statistically uniform manner.  The product coefficients of the four points amount to barycentric coordinates with a sum of one.  The expression

 (1-r3)*P3 + r3*P4

would distribute uniformly along the line between P3 and P4 using r3 = rand.  The quantity

 (1-r2)*P2 + r2*((1-r3)*P3 + r3*P4)

would then distribute uniformly within the triangle P2P3P4 with r2 = rand^(1/2) since the area of a triangle varies as the square of its height.  Finally the full expression with r1 = rand^(1/3) gives a uniform distribution within the tetrahedron P1P2P3P4 since its volume is proportional to the cube of its height.

Roger Stafford