Re: Random points in triangle
- To: mathgroup at smc.vnet.net
- Subject: [mg111644] Re: Random points in triangle
- From: Ray Koopman <koopman at sfu.ca>
- Date: Sun, 8 Aug 2010 07:23:43 -0400 (EDT)
- References: <i3jc3s$cuh$1@smc.vnet.net>
On Aug 7, 3:21 am, Bill Rowe <readn... at sbcglobal.net> wrote:
> On 8/6/10 at 6:55 AM, stev... at ROADRUNNER.COM (S. B. Gray) wrote:
>
>> I was looking for a simple way to place random points inside a
>> triangle with uniform distribution. Here's a good way:
>> newtri := Module[{x},
>> ptri = RandomReal[{-5, +5}, {3, 2}]; tredg = Subsets[ptri, {2}];
>> ]
>> newpts[nump_] := Module[{wts},
>> inpoints = {}; Do [ wts = RandomReal[GammaDistribution[1, 2], 3];
>> wts = wts/Total[wts]; newin = Total[ptri*wts]; inpoints =
>> Append[inpoints, newin], {nump}];
>> ]
>> shotri := Module[{x},
>> Graphics[{Blue, Line[tredg], Red, Point[inpoints]}, ImageSize -> 500]
>> ]
>> The same idea works for points in a tetrahedron; they will be
>> uniformly distributed if you use args such as
>> GammaDistribution[.6,.1].
>
> It seems to me you are making something that should be
> reasonably simple complex and doing so in a manner where you
> aren't truly getting your stated desired result. Here, I take
> your "random points inside a triangle with uniform distribution"
> to mean points uniformly distributed over a specified triangle.
>
> GammaDistribution[1, 2] is the same as ExponentialDistribution[2]
>
> The sum of n independent items drawn from
> ExponentialDistribution[b] will be distributed as
> GammaDistribution[n, b]
>
> The ratio x/(x+y) will have a beta distribution when x and y are
> independent values drawn from gamma distributions.
>
> A uniform distribution can be regarded as a special case of the
> beta distribution with a = b = 1.
>
> So, your variable newwin which is the product of two independent
> things each coming from a beta distribution will certainly not
> be uniformly distributed.
Try this a few times. It looks pretty uniform to me.
p = RandomReal[NormalDistribution[0,1],{3,2}];
r = #.p / Total[#,{2}] & @
RandomReal[ExponentialDistribution[1],{5000,3}];
ListPlot[{p,r},PlotRange->All,AspectRatio->1,Frame->True,Axes->None,
PlotStyle->{{Red,PointSize[.02]},{Black,PointSize[.005]}}]
> However, the actual distribution for
> newwin *might* be close enough to a uniform distribution for
> your purpose. I haven't taken the time to work out the actual
> distribution of newwin.
>
> But rather than going through all of the above and trying to
> work out the exact distribution for newwin, there is a simpler
> approach which is guaranteed to give points uniformly
> distributed over a given shape.
>
> To select random points over some planar shape and ensure the
> points are uniformly distributed over that shape, simply
> generate a pair of uniform deviates using RandomReal[{-1,1}, 2].
> That gives you uniformly distributed points over a square. Now
> all you need do to get uniform points over your desired shape is
> to drop points that lie outside your desired shape.
>
> For simplicity, assume your desire is to have points distributed
> uniformly over a unit circle. Then the values returned by:
>
> Cases[RandomReal[{-1, 1},{1000,2}],_?(Norm@#<=1&)]
>
> will be a set of points uniformly distributed over a circle with
> radius 1
>
> This can be readily see by doing
>
> ListPlot[Cases[RandomReal[{-1, 1}, {1000, 2}], _?(Norm@# <= 1 &)],
> AspectRatio -> 1]
>
> The same idea can be extended for solids. For example,
>
> Cases[RandomReal[{-1, 1},{1000,3}],_?(Norm@#<=1&)]
>
> will be a set of points uniformly distributed through out the
> volume of a sphere with radius 1.
>
> Here, I chose a sphere and circle as examples because it is very
> simple to determine whether a randomly selected point lies
> inside or not. But the idea of doing this selection is valid no
> matter what shape is desired.
>
> Also, by using Cases, I am not controlling how many points I get
> inside the desired area/volume. If you need to have a specific
> number of points inside your desired shape, you could do
> something like:
>
> In[12]:= Table[pt = RandomReal[1, 2];
> While[Norm[pt] > 1, pt = RandomReal[1, 2]]; pt, {5}]
>
> Out[12]= {{0.337431, 0.0873031}, {0.25394, 0.927761}, {0.566057,
> 0.391587}, {0.497942, 0.441549}, {0.0179655, 0.854459}}
>
> which gives you a list of 5 points inside or on the circle with
> radius 1.
>
> And of course it should go with out saying, larger or smaller
> shapes can easily be obtained by multiplying your selected
> points by the appropriate scale factor.
>
> There is one other aspect of this approach that should be
> mentioned. There are more random numbers being generated than
> are actually used. If the desired shape is a very small fraction
> of unit square, there will be a large number of the generated
> points that are simply thrown away, causing a loss of
> efficiency. That is, doing
>
> Cases[RandomReal[{-10, 10},{100000,3}],_?(Norm@#<=1&)]
>
> will generate approximately the same number of points in a
> circle with radius 1 as
>
> Cases[RandomReal[{-1, 1},{1000,3}],_?(Norm@#<=1&)]
>
> at the cost of doing about 100 times as much work.
>
> As long as the area/volume of the desired shape is reasonably
> close to a square/cube then this shouldn't be much of a problem.
> In fact, this method may out perform your approach since it is
> takes less computation to generate a set of uniform deviates
> that it does to generated deviates from other distributions.