[SciPy-User] random points within an ellipse
George Nurser
gnurser@gmail....
Tue Aug 10 11:03:15 CDT 2010
Hi,
Perhaps
#!/usr/bin/env python
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)
ax.scatter(xv,yv)
ax.axis('equal')
plt.show() # Figure attached
HTH. George.
On 10 August 2010 16:29, Sebastian Haase <seb.haase@gmail.com> wrote:
> 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@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@scipy.org
>> http://mail.scipy.org/mailman/listinfo/scipy-user
>>
>>
> _______________________________________________
> SciPy-User mailing list
> SciPy-User@scipy.org
> http://mail.scipy.org/mailman/listinfo/scipy-user
>
More information about the SciPy-User
mailing list