"Bruno Luong" <firstname.lastname@example.org> wrote in message <email@example.com>... > Really cool Roger! I would love if you could share the justification of this formula.
"james bejon" wrote in message <firstname.lastname@example.org>... > 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:
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.