Estimating π using Monte Carlo simulationsTue 05 January 2021
The Monte Carlo method
Monte Carlo methods are a group of computational algorithms that rely on repeated random sampling to obtain numerical results. The underlying concept is to use randomness to solve problems that might be deterministic in principle. They are often used in physical and mathematical problems and are most useful when it is difficult or impossible to use other approaches. Monte Carlo methods are mainly used in optimization problems, numerical integration, and to generate draws from a probability distribution.1
Estimating the value of \(\pi\)
Consider the circle in the diagram. If we generate a sufficiently large number of \((x, y)\) points in the domain \(-1 \leqslant x \leqslant 1\) and \(-1 \leqslant y \leqslant 1\), then the fraction of points that fall inside the circle will be equal to the ratio between the area of the circle and area of the square.
A point falls inside the circle if it complies to the following condition:
By counting the number of random points that fall within the circle we can estimate \(\pi\) as:
We can thus devise an algorithm that will:
- Initialize a counter of the points that fall inside the circle (
- Generate a random point \((x, y)\).
- Check if the point falls within the circle by evaluating \(x^2 + y^2 \leqslant R^2\).
- If the abose condition is true then increment the counter
- Calculate \(\pi\) .
- Repeat from step 2 for the desired number of runs.
In principle, the higher the number of runs the better the estimation will be (closer to the real value of \(\pi\)).
Putting the above algorithm to pratice we get the following code2:
import random def run_simulation(runs=1000): simulation =  circle_points = 0 for i in range(runs): # Create a (x, y) point at random x, y = random.uniform(-1, 1), random.uniform(-1, 1) # Check if the point falls within the circle if x**2 + y**2 <= 1: circle_points += 1 # Calculate current pi value pi_estimate = 4 * circle_points/(i+1) # Save simulation step simulation.append([x, y, pi_estimate]) # Print final pi value print('Result: Pi =', simulation[-1]) return simulation sim = run_simulation(runs=10000)
Here is the result of the simulation from 1000 runs to 100000 runs. At 100000 runs the estimation (\(\pi = 3.142720\)) is already very close the actual value (\(\pi = 3.141593\)).
And if you are curious, here is how the estimation approaches the real value from the first run to 10000th run.
A final remark: We do not need to perform the simulation above using a full circle. By selecting a quarter of a circle and the corresponding square we could achieve the same result with a fourth of the computational effort.