There's an interesting way to use Monte Carlo simulation to find the value of pi. Let's imagine a circular target (without a bullseye this time), with a square that just barely encloses the circular target:
[image in blog post]
As we found in the previous blog post, the probability of hitting this target is:
$p_{hit} = \frac{TargetArea}{TotalArea} = \frac{\pi r_{t}^{2}}{(2r_{t})^{2}} = \frac{\pi}{4}$
When we run the Monte Carlo simulation, we'll choose random numbers between -1 and 1 in both the x and y direction, from a uniform random distribution. The scoring per hit can be reasoned as follows.
If the target is hit, p_hit is 1, which according to the equation above gives:
pi = 4
If the target is not hit, p_hit is 0, so:
pi = 0
So for each target hit, we'll add 4 to our total score, and for each target missed, we'll add 0 to our total score. If we take the average score over each Monte Carlo iteration, this should give us a decent approximation of pi.
The following code implements a Monte Carlo simulation for this problem.
In [13]:
import random
import pandas as pd
from numba import jit
import matplotlib.pyplot as plt
# This is an array, so we can show the improvement
# as increasing number of histories are used
num_hists = [1e1, 1e2, 1e3, 1e4, 1e5, 1e6, 1e7, 1e8, 1e9, 1e10]
# Monte Carlo simulation function. This is defined as
# a function so the numba library can be used to speed
# up execution. Otherwise, this would run much slower.
@jit
def MCHist(n_hist):
score = 0
for n in range(1, n_hist):
x = random.uniform(-1, 1)
y = random.uniform(-1, 1)
# Check if the point falls inside the target
if (x**2 + y**2) <= 1:
# If so, give it a score of 4
score += 4
return score
# Run the simulation, iterating over each number of
# histories in the num_hists array
results = {}
for n in num_hists:
results[n] = MCHist(n) / n
# Show the results in a table
df = pd.DataFrame.from_dict(results, orient="index")
df
Out[13]:
Notice that we can get the first three digits after the decimal place correct using only 1e6 histories, and the first four digits right using 1e8 histories. However, the dropoff in accuracy after adding an order of magnitude may be steep. Even using 1e10 histories only gives us the first 5 digits after the decimal. It's clear that this method of calculating pi is not sustainable past the first few digits.
This demonstrates a rule in Monte Carlo simulation called N^(-1/2). Essentially, getting an extra digit of accuracy in your results requires that you run 100 times more histories than before. Depending on the accuracy required for a problem, the number of histories can get prohibitively large.