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