%pylab inline

nb_pts = 5000
px = rand(nb_pts) # x coordinates
py = rand(nb_pts) # y coordiantes
inside = find((px**2 + py**2) <= 1.0)
nb_inside = len(inside)
ratio = float(nb_inside)/nb_pts
print("pi estimate: %5.3f" % (ratio * 4))

axis('scaled')
plot(px, py, '.')
plot(px[inside], py[inside], 'r.')

# quarter of a circle
cx = linspace(0,1, 100)
cy = sqrt(1 - cx**2)
plot(cx, cy, 'k-', linewidth=2)

import sys, IPython, numpy, matplotlib
print 'Python %d.%d.%d' % (sys.version_info.major, sys.version_info.minor, sys.version_info.micro)
print 'Ipython', IPython.__version__
print 'NumPy', numpy.version.version
print 'MPL', matplotlib.__version__