monte carlo 101

monte_carlo_101.py

import numpy as np
from matplotlib import pyplot as plt

# helper function to check if (x, y) is in red section
# x and y may be a single point or a numpy array
def is_in_quarter_section(x, y):
    return x*x + y*y <= 1

def approximate_pi(sample_size=10000):
    # initialization
    random_seed_initializer = 1
    low = 0
    high = 1
    np.random.seed(random_seed_initializer)
    xs = np.random.uniform(low, high, (sample_size,))
    ys = np.random.uniform(low, high, (sample_size,))

    # pi approximation
    approxi_pi = is_in_quarter_section(xs, ys).sum()/sample_size * 4
    print(f'PI is about {approxi_pi}')

    # draw
    point_colors = ['r' if is_in_quarter_section(zip_xy[0], zip_xy[1]) else 'b' for zip_xy in zip(xs, ys)]
    plt.scatter(xs, ys, c=point_colors)
    plt.show()