[SciPy-User] random points within an ellipse
Tue Aug 10 11:03:15 CDT 2010
import numpy as np
from numpy.random import random_sample as random
a = 1.; b = 0.5
num = 1e4
x = -a + 2*a*random(num)
y = -b + 2*b*random(num)
ra2 = 1./(a*a)
rb2 = 1./(b*b)
# create boolean mask
inside = x*x*ra2 + y*y*rb2 <1.
xv = x[inside]
yv = y[inside]
import matplotlib.pyplot as plt
fig = plt.figure()
ax = fig.add_subplot(111)
plt.show() # Figure attached
On 10 August 2010 16:29, Sebastian Haase <firstname.lastname@example.org> wrote:
> 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.
> Sebastian Haase
> On Sun, Aug 8, 2010 at 11:03 PM, nicky van foreest <email@example.com> 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
> SciPy-User mailing list
More information about the SciPy-User