Day 9: Monte carlo - π

I guess I’ll introduce more the one randomised algorithm in this series, and that’s not only because I love probabilistic approach. Randomised simulation is often the best/only way to solve otherwise intractable problems.

if you are mathematician, excuse my sloppy wording; randomised simulation should be understood as a random sampling from space of simulations under given distribution

How can we estimate π if the only tool we have at disposal is a good random number generator? When we choose a random coordinate (x, y) in range (-1, 1) and each point has equal chance to be chosen, the probability to hit a circle with unit radius is

Having sufficiently large set of points [and a good generator] we can get as close as we want according to Chebyshev’s inequality.


def pi(n, batch=1000):
t = 0
for i in range(n // batch):
p = np.random.rand(batch, 2)
p = (p * p).sum(axis=1)
t += (p <= 1).sum()
return 4 * t / n


> pi(10 ** 8)
One clap, two clap, three clap, forty?

By clapping more or less, you can signal to us which stories really stand out.