In this post I will look at a new variant of an old way of approximating pi using random numbers, a so-called “Monte Carlo” method. Yes, Monte Carlo as in the casino.
There is another way of doing this simulation using darts => link and even more different method called “Count Buffons needle” => link. In this post I wanted to change it up a bit, because others have already used the darts analogy and it is nice to do something in a little different way.
Background
I was at my summer-house this week and thought that using the idea of an apple tree would make the core idea of the Monte Carlo method even more clear and we can all eat apple pi(e) afterwards. By the way, isn’t special in this way, we could have approximate using the same method. Maybe even more interesting, we can approximate integrals in 2 or higher dimensions, see link.
The classic way of calculating using random numbers is to throw darts completely randomly within a square. Inside this square is the board (also called a “clock”) with its radius half of the square’s side. We need to throw many darts completely at random. The ones sticking to the disk are counted and the number of darts sticking to the square are also counted. By using the ratio between the two areas, we can calculate . The details how to do this are described below.
Make like a tree and get outta here
Question: Why not show the idea using the darts idea, isn’t darts good enough for you? Well, I wanted something different and I thought that the idea of throwing darts is not random. A person throwing darts will never achieve a completely random distribution, because that is not the goal of darts, so using darts to explain the idea is not perfect because uniformly random distribution of points is the most important idea in this simulation. Instead, using apples that are dropped from a height is almost as random as it gets, it feels more like Monte Carlo.
Besides, a apple tree is very common, we associate it with Newton and apple pie. We can also instead of counting the number of apples, we can actually weigh them. The answer will be the same and there are usually lots of them, as opposed to darts.
Food for thought
We begin by placing a square plastic sheet under the tree, so the sheet is just inscribed in the shadow of the tree (if the sun is above the tree). We ignore that the tree has a stem, the calculations will be different, but by argument we can see that the result will be the same.
Imagine we gently shake the tree, the apples will land on the sheet. When we are finished, the sheet is weighed and recorded and the other apples are also weighed. That’s it, this is all we need to do. How is this approximating ?
The math
We define the radius of the disc and the side of the square with r and d respectively. Using Pythagoras’ theorem, we can get the side of the square, d, expressed as r. For simplicity, say d = 2k. We then have , so this means that . We need this to calculate the area of the square.
We define the area of the shadow of the tree as and the area of the square as . We have and . = which is simply . = .
So if we let r=1 then . But, we don’t know the areas, how can we get if we don’t know the left hand side?
– The ratio of the areas are the same as the ratio of the number of apples in the areas. If the area is twice as large, then we will have twice as many apples. The problem is that we need many apples. Millions of them. This is why I use a Matlab script to do it for me.
The coding
We need to randomly place points on a disc. In the dart version we generate random points in a square and find how many points are within a inscribed circle, this is trivial. In this version it is the opposite, which presents an interesting problem.
To generate the random points we can simply use polar coordinates. The problem with this is that the number of points are not evenly distributed over the disk. To convince ourselves, we can plot the density by counting the number of points within a radius r for a couple of points.
The reason the points are denser in the center of the disk is a known problem and the solution is to map the radius to => link to mathworld.
r = 2*rand(n, 1); r = sqrt(r); %mapping to uniform distributed samples theta = 2*pi*rand(n, 1); x = r.*cos(theta); y = r.*sin(theta); %pick out the points on the square mask = abs(x) < 1 & abs(y) < 1; xs = x(mask); ys = y(mask);
The variable mask is a binary vector, where the ones are the points within the square and the zeroes are the points outside. By using the number of points to approximate the area, we can use the equation of which we derived above to get our approximation. is thus 2*length(x)/length(xs).
An interesting part of my version is that we can use physical apples and weigh them instead of counting them, if we would like. The weight doesn’t add anything to the ratio. In fact, if you try this at home, you just need a couple of thousand apples to get a good approximation of .
As we can see, the approximation is converging very slowly. It takes about a million apples to get only three decimals. This is about the same as the dart version.
How do you like those apples?
We used a new variant of an old Monte Carlo simulation to approximate . We saw that the problem description, although the same as for the dart version provide us with a additional task of mapping random variable in polar coordinates and it also gives us a nicer connection to the real world. Much like the dart version it takes many samples to get a good approximation. We also discussed that it is possible to approximate other constants such as or even integrate n-dimensional functions using the Monte Carlo method.
Illustrated, conceived and written by Matz JB 2014