[SciPy-User] random points within an ellipse
Tue Aug 10 10:29:34 CDT 2010
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.
On Sun, Aug 8, 2010 at 11:03 PM, nicky van foreest <firstname.lastname@example.org> wrote:
> 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.
> #!/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):
> 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
More information about the SciPy-User