How to estimate the value of π…(part 2)

In my last blog entry, I had ended by introducing the idea why a Monte-Carlo simulation that estimates a \pi without estimating the error involved in the simulation gives an incomplete picture of the simulation. The Monte Carlo Standard Error (MCSE) is a measure of the accuracy of the simulation. The lower the MCSE, the more accurate the simulation.

Let \widehat{\pi} be the Monte Carlo-estimated value of pi and let \widehat{\sigma} be the corresponding standard deviation of the simulation. The MCSE of the simulation would be \frac{\widehat{\sigma_M}}{n} where n would be the size of the generated scenarios (I belabor this point but we generated 2n numbers so we could get n-pairs of numbers and thus we have n scenarios)

So, let us tweak the previous program to change the function estimate_pi()

1def estimate_pi(n=1000000):
2    nums = 2 * np.random.random_sample((n, 2)) - 1.0
4    insides = np.array([1 if inside_circle(x, y) else 0 for x, y in nums])
5    mean = np.mean(insides)
6    sd = np.std(insides)
7    return 4 * mean, 4 * sd/sqrt(n)

This generates both the estimate of \pi and the MCSE. For example, I ran the program with n=1000000, and n=4000000, which gave the estimates and MCSE values of \widehat{\pi}=3.139252 and \widehat{\sigma_M}=0.0016438080424721131 for n=1000000 and \widehat{\pi}=3.141344 and \widehat{\sigma_M}=0.0008211780978667174 for n=4000000