[SciPy-User] random points within an ellipse

[Sebastian Haase]

Hi,
in case you need speed...
if you generate (and test and reject) always 1000 xy points at a time,
you will be MUCH (100x ?) faster.
Try using vector operations where possible.

Regards,
Sebastian Haase


On Sun, Aug 8, 2010 at 11:03 PM, nicky van foreest <vanforeest at gmail.com> wrote:
> Hi,
>
> Good idea to include the graphical result.
>
> I must admit that the results of the latest mail by David do not
> satisfy my "random intuition". There is a thick cloud around the
> origin, which seems unnatural to me.
>
> Just for fun I implemented a rejection method that seems to give good
> results. From an algorithmic perspective my implementation is far from
> optimal (BTW, I welcome feedback on how to speed up code like this.)
> However, the rejection method does not seem that bad, roughly 20% of
> the numbers is rejected, and there is no need to compute cosines and
> the like. The resulting graph seems ok to me.
>
> I modified the example of DG somewhat for the example below.
>
> Nicky
>
> #!/usr/bin/env python
>
> import numpy as np
> from numpy.random import random_sample as random
>
> a = 1.; b = 0.5
>
> def inEllipse(x,y):
>    return (x/a)**2 + (y/b)**2 < 1.  # taking sqrt is unnecessary
>
> num = 1e4
> x = -a + 2*a*random(num)
> y = -b + 2*b*random(num)
>
> sample = []
> reject = 0
> for xx, yy in zip(x,y):
>    if inEllipse(xx,yy):
>        sample.append([xx,yy])
>    else:
>        reject += 1
>
> print reject
> sample = np.array(sample)
>
> import matplotlib.pyplot as plt
>
> fig = plt.figure()
>
> ax = fig.add_subplot(111)
>
> ax.scatter(sample[:,0], sample[:,1])
>
> plt.show() # Figure attached
>
> _______________________________________________
> SciPy-User mailing list
> SciPy-User at scipy.org
> http://mail.scipy.org/mailman/listinfo/scipy-user
>
>