A Monte Carlo simulation that estimates pi. Just like on Wikipedia! http://en.wikipedia.org/wiki/Monte_carlo_method

pi_montecarlo.py

import numpy as np
import matplotlib.pyplot as plt
from matplotlib import animation


def new_figure():
    figure = plt.figure(122, figsize=(10.0, 5.0))
    figure.add_subplot(121)
    plt.axis('equal')

    xs = np.arange(0, 1, .01)
    ys = np.sqrt(1 - xs ** 2)

    plt.plot(xs, ys, lw=5.0, c='#444444')
    plt.fill_between(xs, ys, color='#eeeeee')
    annotations = ['Darts:', 'Hits:', 'Ratio:', '$\pi$:']
    for i, a in enumerate(annotations):
        if a == '$\pi$:':
            plt.text(1.2, 1.0 - 0.1 * i, a, size=24.0)
        else:
            plt.text(1.2, 1.0 - 0.1 * i, a, size=18.0)
    for i in range(4):
        plt.text(1.2 + .5, 1.0 - 0.1 * i, 0, size=18.0, color='k')
    return figure


def throw_dart(point, counter):
    x, y = point
    if np.sqrt(x ** 2 + y ** 2) <= 1.0:
        plt.plot(x, y, marker='o', c='#cc0000', alpha=0.5)
        counter['hits'] += 1.0
    else:
        plt.plot(x, y, marker='o', c='#1100cc', alpha=0.5)
        counter['misses'] += 1.0


def clear_stats():
    print len(plt.gca().texts)
    if len(plt.gca().texts) > 5:
        del plt.gca().texts[-6:]


def update(counter):
    hits = counter['hits']
    misses = counter['misses']
    total = hits + misses
    if misses != 0:
        ratio = hits / total
    else:
        ratio = 0
    pi = 4.0 * ratio
    new_stats = [int(total), int(hits), ratio, pi]
    for i, a in zip(range(4), new_stats):
        j = -4 + i
        plt.gca().texts[j].set_text(a)


def get_point(n):
    i = 0
    while i < n:
        yield np.random.random((1, 2))
        i += 1


def animate(_, points, counter):
    point = points.next()
    throw_dart(point[0], counter)
    update(counter)


def main(*args):
    figure = new_figure()
    counter = {'hits': 0.0,
               'misses': 0.0}
    points = get_point(100000)

    if 'animate' in args:

        anim = animation.FuncAnimation(figure,
                                       animate,
                                       frames=10000,
                                       interval=1,
                                       fargs=(points, counter),
                                       blit=False)
    else:
        for point in points:
            throw_dart(point[0], counter)
        update(counter)


    plt.show()

if __name__ == '__main__':
    main()