[SciPy-User] random points within an ellipse

[David Goldsmith]

Ben:

Nicky brings up a good point: you need to be clear on precisely how you want
the points to be distributed; I was thinking more about my final "working"
example in which the area density appears to be inversely related to r
(which I attributed to the size of the points being more significant with
decreasing r) and I think (I'm too lazy/have too many other priorities to
prove) that the true explanation is that my distribution is uniform with
respect to r (and axially symmetric, which you almost certainly want), which
makes the expected number of points per unit area decrease as 1/r (again,
this is a conjecture, I'm not claiming to have proven it).  If you want a
distribution in which the expected number of points per unit area is
independent of r (and theta), i.e., a "unit-area-uniform" distribution, I
think either Alan's or Erin's methods are what you want.

DG

On Thu, Aug 5, 2010 at 2:16 PM, nicky van foreest <vanforeest at gmail.com>wrote:

> Hi,
>
> I actually wonder whether this is the type of distribution you want to
> see. Lets apply the method described below to a circle (that is just a
> special type of ellipse) with unit radius and center at (0,0). With p
> =1/2 you will hit the segments of the circle  at the left of the line
> such that x=-1/2 or at the right of the line with x =1/2. However the
> segments above the line with y=1/2 or below the line with y = -1/2
> will be hit with a smaller probability. Thus, the distribution of
> points on the circle produced by the method below is certainly not
> rotationally invariant.
>
> I recall having read three ways to produce a random point in a circle,
> but I forgot where (perhaps it can be found in Feller, or the book by
> Jaynes), and each of the three methods gives different distributions
> with different properties.
>
> Nicky
>
> On 5 August 2010 19:53, Ian Nimmo-Smith <iannimmosmith at gmail.com> wrote:
> > Ben
> >
> > You could try the following to generate a uniform distribution in the
> > interior of an ellipse with semi-axes of length (a,b). It scales the
> > x-coordinate of a random point (x,y) in the circular disc of radius 1 by
> a.
> > This will have the correct marginal density of the x-coordinate for the
> > desired ellipse. Conditional on x, the desired y distribution is uniform
> on
> > the chord of length 2*b*sqrt(1-(x/a)**2).  (It's a bit wasteful in that
> it
> > uses 3 calls to random() but discards the y-coordinate of the random
> point
> > in the disc. I don't think you can utilise that to generate a second
> point
> > in the ellipse because y and x are not independent.)
> >
> > Ian
> > ---
> >
> > import sympy
> > import numpy as np
> > import scipy as sc
> > from numpy.random import random_sample as random
> >
> > def random_uniform_in_disc():
> >     # returns a tuple which is uniform in the disc
> >     theta = 2*np.pi*random()
> >     r2 = random()
> >     r = np.sqrt(r2)
> >     return np.array((r*np.sin(theta),r*np.cos(theta)))
> >
> > def random_uniform_in_ellipse(a=1,b=1):
> >     x = a*random_uniform_in_disc()[0]
> >     y = b*np.sqrt(1-(x/a)**2)*(1-2*random())
> >     return np.array((x,y))
> >
> > import matplotlib.pyplot as plt
> > fig = plt.figure()
> > ax = fig.add_subplot(111)
> > sample = np.array([random_uniform_in_ellipse(a=2,b=1) for i in
> > np.arange(10000)])
> > ax.scatter(*sample.T)
> > plt.show()
> >
> >
> > On 5 August 2010 16:06, Benjamin Root <ben.root at ou.edu> wrote:
> >>
> >> On Thu, Aug 5, 2010 at 1:54 AM, David Goldsmith <
> d.l.goldsmith at gmail.com>
> >> wrote:
> >>>
> >>> Whoops, sorry, in the example, the call to rand_in_ellipse should have
> >>> 4./3 as a single argument, i.e., the line should look like
> >>>
> >>> sample = np.array([rand_in_ellipse(4./3) for i in np.arange(10000)])
> >>>
> >>>
> >>>
> >>>
> >>>
> >>> DG
> >>>
> >>> On Wed, Aug 4, 2010 at 11:50 PM, David Goldsmith
> >>> <d.l.goldsmith at gmail.com> wrote:
> >>>>
> >>>> For the ellipse:
> >>>>
> >>>> import numpy as np
> >>>> from numpy.random import random_sample as random
> >>>>
> >>>> pi = np.pi
> >>>> def rand_in_ellipse(a, b=1, offset=0):
> >>>>     angle = 2*pi*random() - offset
> >>>>     x = a * random() * np.cos(angle)
> >>>>     y = b * random() * np.sin(angle)
> >>>>     return np.array((x,y))
> >>>>
> >>>> import matplotlib.pyplot as plt
> >>>>
> >>>> fig = plt.figure()
> >>>>
> >>>>
> >>>>
> >>>>
> >>>>
> >>>> ax = fig.add_subplot(111)
> >>>>
> >>>>
> >>>> sample = np.array([rand_in_ellipse() for i in np.arange(10000)])
> >>>> ax.scatter(sample[:,0], sample[:,1])
> >>>>
> >>>>
> >>>>
> >>>>
> >>>>
> >>>> plt.show() # Figure attached
> >>>>
> >>>> DG
> >>
> >> Unfortunately, that does not produce a uniform sample.   It is heavily
> >> biased along the major and minor axes.
> >>
> >> Ben Root
> >>
> >> _______________________________________________
> >> SciPy-User mailing list
> >> SciPy-User at scipy.org
> >> http://mail.scipy.org/mailman/listinfo/scipy-user
> >>
> >
> >
> > _______________________________________________
> > SciPy-User mailing list
> > SciPy-User at scipy.org
> > http://mail.scipy.org/mailman/listinfo/scipy-user
> >
> >
> _______________________________________________
> SciPy-User mailing list
> SciPy-User at scipy.org
> http://mail.scipy.org/mailman/listinfo/scipy-user
>



-- 
Mathematician: noun, someone who disavows certainty when their uncertainty
set is non-empty, even if that set has measure zero.

Hope: noun, that delusive spirit which escaped Pandora's jar and, with her
lies, prevents mankind from committing a general suicide.  (As interpreted
by Robert Graves)
-------------- next part --------------
An HTML attachment was scrubbed...
URL: <http://mail.scipy.org/pipermail/scipy-user/attachments/20100805/90c4e1cc/attachment.html>