Accurate and fast PI in Python

In my last article I show how to use numpy and numba to speed up a naive Monte Carlo method to compute PI. We managed to get > 20x speed up. However the accuracy of PI even with 10 million samples was only accurate 2 decimal places. Rather disappointing :)

In this article I will show you how scientific programmers or numerical people do MC. Using the following we compute PI as 3.1415348.

The key is to use RandomState and two separate seeds to get two minimally correlated sequences of floats to avoid 'x' and 'y' be correlated. Even though numpy is using a very good Mersenne Twister, the previous method use two consecutive random numbers from the random generator. That is not a good idea.

We could have use a single random generator if we generate all the 'x' first followed by 'y'. But that is less flexible. Imagine you want to compute until certain acceptance criteria is meet. In that case you do not know how many 'x's you will need before you start the simulation.