MathGroup Archive 2010

[Date Index] [Thread Index] [Author Index]

Search the Archive

Re: A distribution problem using Mathematica

  • To: mathgroup at smc.vnet.net
  • Subject: [mg111502] Re: A distribution problem using Mathematica
  • From: Daniel Lichtblau <danl at wolfram.com>
  • Date: Tue, 3 Aug 2010 06:40:48 -0400 (EDT)

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


  • Prev by Date: Re: assuming certain properties about variables
  • Next by Date: Re: ReplaceAll and ReplaceRepeated Strange Behavior
  • Previous by thread: Re: A distribution problem using Mathematica
  • Next by thread: Re: A distribution problem using Mathematica