Integrating a semicircle re-visited

Use Monte Carlo integration (with 100 random points) to calculate the value of the integral:

$ I = \int_{-1}^1\sqrt{1-x^2}\mathrm{d}x $

How does this compare to exact answer? (Hint: the integrand is a semicircle of radius 1)

How can you improve the accuracy of your estimate?

Increase the number of points until you get an accuracy comparable (same order of magnitude) as the Riemann sum method with 100 points (which was implemented in the Numerical integration - quick test).

Use the %%timeit notebook magic to compare the calculation times for the Riemann sum method and Monte Carlo method.

Which is more efficient?

Show answer

We want to calculate the area of a semicircle with radius 1. We can adapt the approach used in the Monte Carlo tutorial but, in this case, we use the fact that $P_i = \frac{A_c}{A_r}$ where $A_r$ is a rectangle of length 2 (as the semicircle goes from $-1$ to $1$) and height 1.

import random
import math
    
def estimate_semicircle_area(num_points):
    
    points = []
    hits = 0
    for i in range(num_points):
        # random.uniform(a,b) returns a random number drawn from a uniform distribution from a to b
        x, y = random.uniform(-1,1), random.uniform(0,1)
        # we test if the point is within the circle (using the equation for a circle, X^2+y^2=r^2)
        if x*x + y*y < 1.0:
            hits += 1
    
    probability = hits / num_points
    rectangle_area = 2
    return probability*rectangle_area
    
estimate_semicircle_area(100)

1.58

Note that your estimate will be different as you be using a different set of random numbers.
The exact answer is $\frac{\pi}{2}$. The error on our calculation is

math.pi/2 - estimate_semicircle_area(100)

-0.04920367320510355

To improve the accuracy we can use a larger number of random points:

math.pi/2 - estimate_semicircle_area(1000)

0.04679632679489654

math.pi/2 - estimate_semicircle_area(10000)

-0.0048036732051033315

Increasing the number of points to 10,000 gives an error comparable to the Riemann sum method with 100 integration slices (where the error is 0.002). Let’s use the %%timeit magic to time how long each takes to run

%%timeit
estimate_semicircle_area(10000)

6.23 ms ± 370 µs per loop (mean ± std. dev. of 7 runs, 100 loops each)

%%timeit
rectangular_slice_integral(semicircle, -1, 1, 100)

31 µs ± 1.09 µs per loop (mean ± std. dev. of 7 runs, 10000 loops each)

As each method gives roughly the same accuracy, but the estimate_semicircle_area is 100x smaller, we can deduce that the Monte Carlo method implemented in estimate_semicircle_area is considerably less efficient than the Riemann summation method implemented in rectangular_slice_integral. However the Monte Carlo method is useful for badly behaving systems, as we will see in the lesson exercises.