|
|
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(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
|
|
