
Re: Generate random numbers with fixed sum and different constraints
Posted:
Nov 20, 2012 2:54 PM


"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(m1,n),1./(m1: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 m1 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)*(1r1) ; (1*r1)*(1r2) ; (1*r1*r2)*(1r3) ; (1*r1*r2*r3)*(10) ];
which with the definitions of the points is the same as
X = (1r1)*P1 + r1*((1r2)*P2 + r2*((1r3)*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
(1r3)*P3 + r3*P4
would distribute uniformly along the line between P3 and P4 using r3 = rand. The quantity
(1r2)*P2 + r2*((1r3)*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

