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