Compute pi by MCMC

mc_pi.py

from random import random

"""
Find pi by the Monte-Carlo method.

area of a circle = pi r^2
area of a square = (2r)^2 = 4 r^2

Perform random uniform sampling between -1 and 1.
The proportion of points in the unit circle is:

p = (pi r^2) / (4r^2)

So pi = 4.0 * p
"""

N = 1000000

n_circle = 0
n_total = 0

for i in xrange(N):
    x = random() * 2 - 1
    y = random() * 2 - 1

    # equation of the unit circle
    if x ** 2 + y ** 2 < 1:
        n_circle += 1

    n_total += 1

print "pi = ", 4.0 * n_circle / n_total