[SciPy-User] random points within an ellipse

[Ian Nimmo-Smith]

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
>
>
-------------- next part --------------
An HTML attachment was scrubbed...
URL: <http://mail.scipy.org/pipermail/scipy-user/attachments/20100805/c9964037/attachment.html>