01 August 2007

another way to estimate π

Here's the slowest, most inefficient method known for computing p: Choose an arbitrary initial value x (real) and iterate the mapping x -> (2-x)/(3-4x) for a total of M times. Count the number of iterates that fall within a small interval around zero. For example, let N denote the number of iterates that fall between -u/2 and +u/2 for some small u. Then the limiting value of the ratio uM/N is [π] (for sufficiently large M as u goes to zero).

I'm reminded of Buffon's needle, another experimental method for computing π. If we have a floor with parallel lines drawn on it at regular intervals, and we toss a needle on the floor at random, the average number of lines on the floor which are crossed by the needle is related to π. In particular, it's 2l/(tπ) where l is the length of the needle and t is the spacing of the lines. (For those who have seen this problem before, there's often a restriction that the needle has to be shorter than the distance between the lines; then the average number of lines crossed is just the probability that a line is crossed on any given toss of the needle.) It's possible to solve this for π. It's also possible to cheat and say you did the experiment and got a very good approximation to π, by taking into account the well-known approximation 355/113 for π and setting up the probability appropriately. This is described at the Wikipedia article. I'm not sure why anyone would bother to do this.

Another way to approximate π is to pick points uniformly at random in a square of side 2, in which is inscribed a circle of radius 1; the proportion of the points which are in the circle is π/4. I doubt that anybody can pick points uniformly at random from a square with sufficient accuracy for this to work. A computer isn't bad at picking points at random, though, and this leads to the following method of approximating π, which I've heard of before -- pick X and Y uniformly at random, and independently, from [-1, 1]. (This corresponds to picking a point in the square.) If X2+Y2 ≤ 1 (and thus the point is in the circle) increment a counter. Doing twenty such trials ten thousand times I get the following numbers of hits:

which correspond to estimates of π/4 obtained by just placing a decimal point in front of each of these; we really have π/4 = .7853981634... It's not that surprising that these estimates are a bit spread out. The number of hits is given by a binomial distribution with n = 10000, p = π/4; thus its mean is 2500π (about 7854) and its standard deviation (10000 (π/4) (1-π/4))1/2, or about 41.

However, I can try out the "Atlantic City Method" in the comfort of my living room. (I could try the Buffon's needle method here, since I have the right sort of floor, but I'm not that bored.) Picking a number uniformly at random from [0,1], and carrying out the "Atlantic City Method" as described above with u=.01, N=10000, I get the following values for M:

and again these correspond to approximations of 1/π if you put a decimal point in front; not bad! It seems surprising that increasing u should do this -- we're told that we should let u approach zero -- but we're able to reduce the ratio of the standard deviation of the number of hits to the mean, and thus make the approximations better, by increasing u.

However, if we increase u even more this doesn't hold up; for example if u = 1 then M = 3524 or 3525. The random aspect of the problem seems to go away but the constant we get is not the one we were looking for! (And if I didn't know the value of π already, then I wouldn't be able to assess this.)

Now, why should this method work? The map f: x -> (2-x)/(3-4x) is an example of a Möbius transformation, also known as a fractional linear transformation. This insight is important, because Mobius transformations can be viewed as isometries of the Riemann sphere; it's clear that this particular transformation takes real numbers to real numbers. My first instinct was that it had something to do with the magnitude of the rotation represented by this particular Mobius transformation, but that doesn't seem to be anything nice. However, near zero the projection taking the complex plane to the Riemann sphere very nearly doubles lengths. (See the picture at the Wikipedia article on the Riemann sphere.) So the length of interval u around zero is mapped to an arc of interval just slightly smaller than 2u on the "real circle" -- that is, the projection of the real axis, plus a point at infinity, onto the Riemann sphere. The sequence of iterates x, f(x), f(f(x)), ... is what one might call a "pseudo-random" sequence of points (it's certainly not random, but I'll think of it as if it were) and it lies in the arc in question a proportion N/M of the time; however it also lies in the arc a proportion just slightly less than 2u/(2π) of the time if you believe the psuedo-randomness, since the arc has length slightly less than u and the circle is a circle of radius 1. As u approaches 0 the error in this approximation goes to zero. So we get N/M = u/π, or π = uM/N. (Of course, this is a heuristic argument, and certainly not a proof, but I just wanted to convince myself that it made sense.) What's more, the actual identity of the linear transformation appears to be a red herring; I suspect this works for a rather large class of Mobius transformations. (Not all, though; for example, for some choices of f the sequence x, f(x), f(f(x)), ... approaches one of the fixed points of f.) I'm having trouble identifying other Mobius transformations that actually have this property, though, so I'm not totally confident in saying this.

11 comments:

The needle-dropping method (with the length of the needles exactly equal to the spacing between the lines) pops up in Samuel Delany's Neverÿon series at one point. Better still, in an author's note at the end of a later volume in the same series, he prints a letter from a mathematician who points out that the approximation of pi you get this way converges to the correct number much more slowly than was implied by the fictional text. (If I remember right, to get n degrees of precision, you need on the order of 10^(2n) tosses...)

with the needle length equal to the line spacing, the probability of a needle hitting one of the lines is 2/π, and each such hit is independent. The expected number of hits after t tosses is 2t/π, with variance t(2/π)(1-2/π); what's important is that this is propotional to t. Thus the standard deviation of the number of hits after t tosses is proportional to t^(1/2). The standard deviation of the estimate of 2/π thus obtained is this divided by t, or t^(-1/2). So you're right about the rate of convergence. From what I understand this is pretty typical for Monte Carlo methods in general.

The "why" of Buffon's needle is that it's a very simple demo of how Monte Carlo techniques work (and how they don't work). It's more for a demo or proof-of-concept than anything practical.

For many problems, when no analytical soln exists and asymptotics are lacking, one can throw lots of computer time at it with Monte Carlo methods and get some results.

The 1/sqrt(t) for order of magnitude of error of estimate is pretty standard for these types of things, which just comes from the property of a sample mean as estimator, assuming nice properties of what you're sampling from and independent samples.

The "Mobius method" for computing pi will work for any transformation f(x)=(ax+b)/(cx+d) such that (a-d)^2 = -4b(b+c). It can also be made to work for other transformations, and for counting iterates near points other than zero, but then uM/N approaches Q*pi, where Q is the square root of a rational number. The explanation for this is given at the end of the note on Linear Fractional Transformantions on that same web site.