I flipped through Tarasov's book The World Is Built on Probability and an example introducing Monte Carlo methods (p. 63 ff.) caught my attention. The idea is to calculate approximations of \(\pi\) using the ratio between the area of a quarter circle \(( \frac{1}{4} \pi r^2)\) with radius \(r\) and the area of the quarter circle's bounding square \((r^2)\). $$ \frac{\frac{1}{4} \pi r^2}{r^2} = \frac{\pi}{4} $$
To get some actual numbers to calculate with, we randomly^{*} sprinkle points \(P = (x, y)\) all over the square. Given enough random points, we expect the ratio between the number of points within the quarter circle (\(x^2 + y^2 < r^2\)) and the total number of points to get close to the ratio of the areas.
For an integer square number of points, we can alternatively arrange them with equal spacing like a grid and count as above. To get a better estimate, increase the number of points, either by making them more dense, or by enlarging the square.
Another (boring) way of approximating the area ratio is to calculate the area of the quarter circle: slice the square into vertical strips of equal width and calculate where the centerline \(x = x_p\) of each strip (or really any vertical line within) intersects the circle \(x^2 + y^2 = r^2\) in the first quadrant. $$ y_p = \sqrt{r^2 - x_p^2} $$
Calculate this for each strip, multiply by the width of the strip, and add it all up. This is, of course, just a Riemann sum and as the width of the strips approaches zero, the sum approaches the actual area of the quarter circle.
Now divide by the area of the square \(r^2.\)I ran implementations of the above algorithms in Julia. These are the results (the numbers of the random column fluctuate with each run; correct digits in bold):
# Points | Random | Grid | Rectangles |
---|---|---|---|
10^{2} | 3.08 | 3.16 | 3.141936857900007 |
10^{4} | 3.136 | 3.1428 | 3.141592998024645 |
10^{6} | 3.1436 | 3.141676 | 3.141592653934176 |
10^{8} | 3.14136696 | 3.14159388 | 3.141592653592467 |
10^{10} | 3.141580126 | 3.1415926932 | 3.141592653589404 |
As you can see, the Monte Carlo simulation is a terrible way of approximating \(\pi,\) and the other two aren't really impressive either. There are much better methods. But it does work and is quite fun.
^{*} The method of generating random \(x\) and \(y\) coordinates in \([0, 1)\) given in the book, dividing random integers from a range \([0, 10^N)\) by \(10^N\), is somewhat problematic as it introduces small floating point errors. Not all integer fractions obtained this way can be represented with standard floating-point arithmetic (some aren't dyadic rationals), so they will be rounded to their closest representation and our numbers are no longer uniformly distributed.
Tarasov, L. (1988). The world is built on probability. Mir Publishers.