[SciPy-User] random points within an ellipse

[Benjamin Root]

On Thu, Aug 5, 2010 at 12:53 PM, 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()
>
>
Ian,

This seems to work very nicely.  And I am not too worried about the extra
call to random.  I am only looking to create about 50 or some points at any
given run (I am simulating storm tracks and want to be able to define a
region of storm initiation).

Thank you,
Ben Root
-------------- next part --------------
An HTML attachment was scrubbed...
URL: <http://mail.scipy.org/pipermail/scipy-user/attachments/20100805/d8b51dd1/attachment.html>