Ten years ago, I was in the middle of my MSc degree, titled “Distributed High-Fidelity Graphics Using P2P”. My research revolved around ray tracing – a family of techniques used to produce highly photorealistic images by simulating various aspects of Physics, especially the way light interacts with the various surfaces and materials it comes in contact with.

Ray tracing is very complex and computationally intensive, and it would take an entire book to write anything substantial about it. But, at the heart of modern ray tracing is an approach called Monte Carlo simulation, or stochastic simulation. It’s a fancy way of saying: this problem is too complex and I don’t have an easy way to calculate a solution, so I’ll take a number of random samples and approximate it instead.

In this article, I’ll illustrate the concept by calculating the area of a circle. Since we have an easy formula to calculate the area of a circle (πr²), there’s normally no need to do this via Monte Carlo simulation. However, it’s an easy way to explain the concept, which you can then use to solve more complex problems.

Approximating the Area of a Circle in Theory

Let’s assume that we don’t know how to calculate the area of a circle – either because a formula does not exist, or because it is too complex to derive.

However, we do have a way to determine whether a point lies inside a circle. If the point is represented by the coordinates (x, y), then the point lies inside the circle of radius r if (x² + y²) <= r², assuming that the circle’s centre is located at (0, 0).

In this case, we can take a number of random points in space and compare how many fall inside and outside of the circle. You can see a simple visualisation of this above. Imagine we were throwing darts at a dartboard while blindfolded. The more darts we throw, the more we are filling in the circle, and the ratio of darts inside vs outside the dartboard tells us the overall area we covered.

It’s a little counterintuitive to imagine how a random process leads towards a deterministic solution. In fact, it’s not really deterministic, and the actual result changes with every run; but the approximation does work very well, and the more samples you take, the closer you get to the real solution. Let’s test this out in practice.

Approximating the Area of a Circle in Python 3.8

The code to implement the Monte Carlo approximation of the area of a circle is quite straightforward, as you can see below.

import random

def get_random_coord(radius):
    return random.uniform(-radius, radius)

def calculate_area_monte_carlo(radius):

    total_in_circle = 0

    for i in range(0, 10000000):
        (x, y) = get_random_coord(radius), get_random_coord(radius)
        in_circle = x ** 2 + y ** 2 <= radius ** 2
        
        if in_circle:
            total_in_circle += 1

        if i in [9999, 99999, 150000, 999999, 1500000, 9999999]:
            approximate_area = (2 * radius) ** 2 * total_in_circle / i
            print(f"Area ({i} samples): {approximate_area}")

Essentially, all I’m doing is:

• Defining a helper function get_random_coord() to generate a random coordinate, just so that I don’t have to repeat the same code twice.
• Generating an (x, y) coordinate pair within the square containing the circle.
• Determining whether this (x, y) falls within the circle, using the aforementioned formula (x² + y²) <= r², in which case I increment total_in_circle.
• At specific (arbitrary) intervals (numbers of samples), I calculate the approximate area of the circle and print it out to show how the calculated area gets more accurate, the more samples you take.

The approximate area is calculated as a ratio of samples in circle (total_in_circle) vs total samples both inside and outside of the circle (i). However, we have to adjust that by the area of a square, which is 2r², as explained in this Stack Overflow answer.

In reality, we do know a formula to accurately calculate the area of a circle (πr²), so let’s toss that in as well in order to compare how accurate the results are:

import random
import math

def get_random_coord(radius):
    return random.uniform(-radius, radius)

def calculate_area_monte_carlo(radius):

    total_in_circle = 0

    for i in range(0, 10000000):
        (x, y) = get_random_coord(radius), get_random_coord(radius)
        in_circle = x ** 2 + y ** 2 <= radius ** 2
        
        if in_circle:
            total_in_circle += 1

        if i in [9999, 99999, 150000, 999999, 1500000, 9999999]:
            approximate_area = (2 * radius) ** 2 * total_in_circle / i
            print(f"Area ({i} samples): {approximate_area}")

def calculate_area_formula(radius):
    area = math.pi * radius ** 2
    print("Area (reference):", area)

def main():
    radius = 400
    calculate_area_monte_carlo(radius)
    calculate_area_formula(radius)

if __name__ == "__main__":
    main()

We can now run this program, and after a little patience, we get some results:

$ python3 main2.py
Area (9999 samples): 501874.1874187419
Area (99999 samples): 503793.8379383794
Area (150000 samples): 503940.26666666666
Area (999999 samples): 502511.2225112225
Area (1500000 samples): 502637.2266666667
Area (9999999 samples): 502743.410274341
Area (reference): 502654.8245743669

If we run it again, the results are a little different, but the trend is still there:

$ python3 main2.py
Area (9999 samples): 506930.69306930696
Area (99999 samples): 501861.0186101861
Area (150000 samples): 502101.3333333333
Area (999999 samples): 502941.30294130294
Area (1500000 samples): 502856.1066666667
Area (9999999 samples): 502779.44227794424
Area (reference): 502654.8245743669

This is the essence of the Monte Carlo approach: the more random samples you take, the more you’re filling in the problem space. However, this comes at a cost: the more samples you take, the more time it takes to compute. With enough samples, the error eventually becomes small enough that the approximation is indistinguishable from an accurate solution. Let’s visualise this.

Visualising the Area of a Circle Approximation

We can use the PIL library to quickly generate progressive images of the circle coverage. All we need is to add a few lines to our existing code:

import random
import math
from PIL import Image

def get_random_coord(radius):
    return random.uniform(-radius, radius)

def calculate_area_monte_carlo(radius):
    diameter = radius * 2
    image = Image.new("RGB", (diameter, diameter), (255, 255, 255))
    in_colour = (255, 0, 0) # red
    out_colour = (0, 0, 255) # blue

    total_in_circle = 0

    for i in range(0, 10000000):
        (x, y) = get_random_coord(radius), get_random_coord(radius)
        in_circle = x ** 2 + y ** 2 <= radius ** 2
        
        if in_circle:
            total_in_circle += 1
            colour = in_colour
        else:
            colour = out_colour

        image.putpixel((math.floor(x + radius), math.floor(y + radius)), colour)

        if i in [9999, 99999, 150000, 999999, 1500000, 9999999]:
            approximate_area = (2 * radius) ** 2 * total_in_circle / i
            print(f"Area ({i} samples): {approximate_area}")
            image.save(f'circle_{i}.png')

def calculate_area_formula(radius):
    area = math.pi * radius ** 2
    print("Area (reference):", area)

def main():
    radius = 400
    calculate_area_monte_carlo(radius)
    calculate_area_formula(radius)

if __name__ == "__main__":
    main()

What we’re doing here is:

• Creating an image at the beginning.
• For each sample, we colour the corresponding pixel either red (if it’s in the circle) or blue (if it’s outside the circle). Note that we have to add the radius to x and y because in the image, (0, 0) is not the centre of the circle, but the upper-left corner of the image.
• Saving an image whenever we print an approximate area.

Earlier, we saw how the numbers seemed to be converging towards the accurate solution. If we now run the updated code above, we can verify the progressive coverage of the circle from actual images (click to enlarge):

Applications

As I said earlier, the Monte Carlo approach is not necessary for problems where we have a simple formula, like calculating the area of a circle. However, not all problems have a simple formula.

If you’ve studied calculus, you’ll know that you can integrate a function to calculate the area under its graph. You can do that for a simple function like 5x², but there are more complex functions for which it is impossible to derive an analytical solution. Even if you never studied calculus, you can imagine how hard it must be to calculate the area of an irregular shape, such as your favourite politician. Randomly throwing darts at them and counting how many hit might be one way of working that out.

And if you ever work with ray tracing, you’ll eventually realise that ray tracing is really just another integration problem – one where you have to calculate the area (or contribution of light towards a surface) of something very irregular. A Monte Carlo approximation is a very effective – if computationally expensive – way of solving this otherwise intractable problem.