[SciPy-User] random points within an ellipse

[nicky van foreest]

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
>
>