π by rejection

In the 18th century, Georges-Louis Leclerc, Count of Buffon, proposed an experiment to estimate \pi. The experiment (now somewhat famous as it appears in almost all probability textbooks) consists in dropping randomly matches on a floor made of parallel floorboards and using the number of time the matches land across a junction to estimate \pi.

pi-pie

To perform the Count’s experiment, you do not need a lot of math. You only need to test if

\displaystyle x\leqslant\frac{w}{2}\sin\alpha

with x\sim\mathcal{U}(0,w) and \alpha\sim\mathcal{U}(0,\frac{\pi}{2}) are both uniform random variables, and w is the width of the floorboards. You may remark that we use \frac{\pi}{2} and that it looks like a circular definition until you think that \frac{\pi}{2} radians is 45°, and you can estimate it using other means. Then you start throwing zillions of virtual matches and count the number of intersections, then use a laboriously obtained probability distribution to estimate \pi.

Lucky for us, there’s a variant that does not require the call to sines, not complicated integrals. Just squaring, counting, and a division. Let’s have a look.

Let’s first start by a refresher. The area of a circle of radius r is

A=\pi r^2

If r is conveniently 1, then we have a unit circle and the area is just, now, \pi.

This suggests a variation on Buffon’s experiment to estimate \pi. Instead of throwing matches, let’s say we rather would throw darts at a board. A square board, 2&time;2 units with a unit circle placed in it. Something like this:

board-00

Now suppose I have a really bad dart player that throws dart uniformly randomly on the square board. If I count the numbers of darts that land in the circle vs the total number of darts, I can get an estimate for \pi. Because the area of the circle is \pi and the area of the square is four, the proportion of darts hitting the board inside the circle is \frac{\pi}{4} (asymptotically).

I’m somewhat lazy and did not want to deal with uniformly random numbers on [-1,1], just [0,1], so I cut the board in four. Let us look at the upper-right corner:

board-01

Now, the quadrant (well, yes, that’s the correct name for a quarter of circle, methinks) as an area of \pi/4, and the square is now a unit square, so the quadrant occupies \pi/4ths of the square (since \frac{\frac{\pi}{4}}{1}=\frac{\pi}{4}). And then we play the dart throwin’ game:

If it lands outside,

board-02

I get 0 points, if it lands inside,

board-03

I get one point.

Now, checking if my uniformly random point (x,y)\sim\mathcal{U}((0,1),(0,1)) is inside the circle is easy: if x^2+y^2\leqslant{1}, then it’s inside (and then we count the infinitesimal border as inside).

The code is:

double pitify(size_t rounds)
 {
  size_t in=0;
  for (size_t i=0;i<rounds;i++)
   {
    float x= std::rand()/(double)RAND_MAX;
    float y= std::rand()/(double)RAND_MAX;
     if (sqr(x)+sqr(y)<=1)
      in++;
   }

  return 4*(in/(double)rounds); // estimate pi
 }

Which is straightforward enough. Now, how efficient is it at estimating \pi? Well, not very. Not very much at all.

To measure the efficiency of the quadrant method (let’s call it that even though I’m quite sure it has an official name), let us look at the estimates produced by varying the number of rounds in powers of ten: 1, 10, 10^2, all the way to 10^9. Since it’s random, we could be unlucky and get a very bad result, so let’s repeat every experiments 1000 times. This allows us to get the box-plots:

box-plots-of-digits

The dashed line is at \pi. At rounds=1, we get either 0 or 4, since one dart is either in or outside the circle (and it breaks gnumeric’s box plot routine) and the estimate is 3.08. Then as the number of rounds increases, error goes down. But not terribly fast. Let’s see how close it gets:

Number of darts,
repeated 1000&time;
Estimation of
\pi
Correct
digits
10^0 3.08 0
10^1 3.1568 1
10^2 3.14616 2
10^3 3.140744 2
10^4 3.1404308 2
10^5 3.1416916 3
10^6 3.141599272 5
10^7 3.1415778452 4
10^8 3.14158891828 4
10^9 3.14159450638 5

So the closest to \pi=3.14159265358979\ldots we got after 10^9\times{}10^3=10^{12} darts (and about 4h CPU) is \approx{}3.1415945, which isn’t great. Which is somewhat counter-intuitive since the variance is reduced by a factor 10 each time.

*
* *

So this method is très mauvaise to estimate \pi. Fortunately for us, there are analytic methods to get good estimates of \pi, like Plouffe’s incredible spigot formula, that can yield the nth digit of \pi.

To be continued

2 Responses to π by rejection

  1. […] a previous episode, we looked at how we could use random sampling to get a good estimate of , and found that we […]

  2. […] ancient Egyptians, but can we do better? In previous posts, we have examined how to approximate by rejection and by approximating the area of a circle on a discrete grid. The rejection method showed very slow […]

Leave a Reply

Fill in your details below or click an icon to log in:

WordPress.com Logo

You are commenting using your WordPress.com account. Log Out / Change )

Twitter picture

You are commenting using your Twitter account. Log Out / Change )

Facebook photo

You are commenting using your Facebook account. Log Out / Change )

Google+ photo

You are commenting using your Google+ account. Log Out / Change )

Connecting to %s

%d bloggers like this: