Prev: Surprising FullSimplify result
Next: ReplaceAll and ReplaceRepeated Strange Behavior
From: S. B. Gray on 2 Aug 2010 07:04
I want 3 random variables v1,v2,v3 with uniform distribution from 0 to
1, but modified or normalized so that their sum is a uniform
distribution from 0 to 1.
These variables are for placing a point at a random place inside a
tetrahedron defined by 4 vertex vectors p1,p2,p3,p4. The internal point
is given by p=v1(p1-p4)+v2(p2-p4)+v3(p3-p4) (barycentric coordinates).
The 0 to 1 constraints assure that the point p will be inside.
The closest I have come is this function giving the triplet
frs={v1,v2,v3}:

mul = RandomReal[{0, 1}, {3}];
mul = mul/Total[mul];
frs = Sqrt[RandomReal[{0, 1}, {3}]]*Power[mul, (3)^-1];

which I tested with the histogram nhist: (The +1 avoids trying to access
the 0th element of the list nhist.)

nhist = Table[0, {1000}];
Do [ mul = RandomReal[{0, 1}, {3}];
mul = mul/Total[mul];
frs = Sqrt[RandomReal[{0, 1}, {3}]]*Power[mul, (3)^-1];
ip = IntegerPart[1000 frs][[2]];
nhist[[ip+1]]++, {10000}
];
Print[frs];
ListPlot[nhist]

This ad-hoc method gives a distribution that covers the range 0-1 but is
too heavy in the region 0.3 to 0.7. This would put too many points near
the middle of the tetrahedron. Something tells me there must be a better
and more elegant solution. Any ideas?

Steve Gray

From: Ray Koopman on 3 Aug 2010 06:38
On Aug 2, 4:04 am, "S. B. Gray" <stev...(a)ROADRUNNER.COM> wrote:
> I want 3 random variables v1,v2,v3 with uniform distribution from 0 to
> 1, but modified or normalized so that their sum is a uniform
> distribution from 0 to 1.
> These variables are for placing a point at a random place inside a
> tetrahedron defined by 4 vertex vectors p1,p2,p3,p4. The internal point
> is given by p=v1(p1-p4)+v2(p2-p4)+v3(p3-p4) (barycentric coordinates).
> The 0 to 1 constraints assure that the point p will be inside.
> The closest I have come is this function giving the triplet
> frs={v1,v2,v3}:
>
> mul = RandomReal[{0, 1}, {3}];
> mul = mul/Total[mul];
> frs = Sqrt[RandomReal[{0, 1}, {3}]]*Power[mul, (3)^-1];
>
> which I tested with the histogram nhist: (The +1 avoids trying to access
> the 0th element of the list nhist.)
>
> nhist = Table[0, {1000}];
> Do [ mul = RandomReal[{0, 1}, {3}];
> mul = mul/Total[mul];
> frs = Sqrt[RandomReal[{0, 1}, {3}]]*Power[mul, (3)^-1];
> ip = IntegerPart[1000 frs][[2]];
> nhist[[ip+1]]++, {10000}
> ];
> Print[frs];
> ListPlot[nhist]
>
> This ad-hoc method gives a distribution that covers the range 0-1 but is
> too heavy in the region 0.3 to 0.7. This would put too many points near
> the middle of the tetrahedron. Something tells me there must be a better
> and more elegant solution. Any ideas?
>
> Steve Gray

This will give n points uniformly distributed inside an m-dimensional
polytope whose vertices are given by the rows of a table p whose
dimensions are {m+1,m}:
Table[ #/Tr@# & @ RandomReal[ExponentialDistribution[1], m+1], {n}].p

From: Ray Koopman on 3 Aug 2010 06:39
On Aug 2, 4:04 am, "S. B. Gray" <stev...(a)ROADRUNNER.COM> wrote:
> I want 3 random variables v1,v2,v3 with uniform distribution from 0 to
> 1, but modified or normalized so that their sum is a uniform
> distribution from 0 to 1.
> These variables are for placing a point at a random place inside a
> tetrahedron defined by 4 vertex vectors p1,p2,p3,p4. The internal point
> is given by p=v1(p1-p4)+v2(p2-p4)+v3(p3-p4) (barycentric coordinates).
> The 0 to 1 constraints assure that the point p will be inside.
> The closest I have come is this function giving the triplet
> frs={v1,v2,v3}:
>
> mul = RandomReal[{0, 1}, {3}];
> mul = mul/Total[mul];
> frs = Sqrt[RandomReal[{0, 1}, {3}]]*Power[mul, (3)^-1];
>
> which I tested with the histogram nhist: (The +1 avoids trying to access
> the 0th element of the list nhist.)
>
> nhist = Table[0, {1000}];
> Do [ mul = RandomReal[{0, 1}, {3}];
> mul = mul/Total[mul];
> frs = Sqrt[RandomReal[{0, 1}, {3}]]*Power[mul, (3)^-1];
> ip = IntegerPart[1000 frs][[2]];
> nhist[[ip+1]]++, {10000}
> ];
> Print[frs];
> ListPlot[nhist]
>
> This ad-hoc method gives a distribution that covers the range 0-1 but is
> too heavy in the region 0.3 to 0.7. This would put too many points near
> the middle of the tetrahedron. Something tells me there must be a better
> and more elegant solution. Any ideas?
>
> Steve Gray

Further to my earlier post, this should be a little faster:

#.p/Total[#,{2}]&@RandomReal[ExponentialDistribution[1],{n,m+1}]

From: Daniel Lichtblau on 3 Aug 2010 06:40
S. B. Gray wrote:
> I want 3 random variables v1,v2,v3 with uniform distribution from 0 to
> 1, but modified or normalized so that their sum is a uniform
> distribution from 0 to 1.
> These variables are for placing a point at a random place inside a
> tetrahedron defined by 4 vertex vectors p1,p2,p3,p4. The internal point
> is given by p=v1(p1-p4)+v2(p2-p4)+v3(p3-p4) (barycentric coordinates).
> The 0 to 1 constraints assure that the point p will be inside.
> The closest I have come is this function giving the triplet
> frs={v1,v2,v3}:
>
> mul = RandomReal[{0, 1}, {3}];
> mul = mul/Total[mul];
> frs = Sqrt[RandomReal[{0, 1}, {3}]]*Power[mul, (3)^-1];
>
> which I tested with the histogram nhist: (The +1 avoids trying to access
> the 0th element of the list nhist.)
>
> nhist = Table[0, {1000}];
> Do [ mul = RandomReal[{0, 1}, {3}];
> mul = mul/Total[mul];
> frs = Sqrt[RandomReal[{0, 1}, {3}]]*Power[mul, (3)^-1];
> ip = IntegerPart[1000 frs][[2]];
> nhist[[ip+1]]++, {10000}
> ];
> Print[frs];
> ListPlot[nhist]
>
> This ad-hoc method gives a distribution that covers the range 0-1 but is
> too heavy in the region 0.3 to 0.7. This would put too many points near
> the middle of the tetrahedron. Something tells me there must be a better
> and more elegant solution. Any ideas?
>
> Steve Gray

You can move along edges from first to second, second to third, and
third to fourth vertices to get to a point. To do so with equal
probability, weight so that first move, for random 0<r1<1, goes a
distance that covers r1 of the volume. This means, I think going
r1^(1/3) of the distance along that edge. Then for 0<r2<1, go distance
that covers r2 of the area of the appropriate cross sectional triangle,
etc. Code below purports to do all this.

randomTetPoint[verts_] := Module[
{dirs = -Apply[Subtract, Partition[verts, 2, 1], 1],
vals = RandomReal[{0, 1}, 3]},
vals = Reverse[Flatten[MapIndexed[#^(1/#2) &, vals]]];
vals = Rest[FoldList[Times, 1, vals]];
First[verts] + vals.dirs
]

Quick visual sanity test:

vertices = {{0, 0, 0}, {1, 0, 0}, {0, 1, 0}, {0, 0, 1}};
rndpts = Table[randomTetPoint[vertices], {2000}];
ListPointPlot3D[rndpts, BoxRatios -> {1, 1, 1}]

Daniel Lichtblau
Wolfram Research

 | 
Pages: 1
Prev: Surprising FullSimplify result
Next: ReplaceAll and ReplaceRepeated Strange Behavior