Calculating pi with a random number generator

It's π day (3/14), so let's estimate π, the ratio of circle's diameter to its circumference, with a random number generator. This post is a blatant ripoff of Rhett Allain's post on the occasion of π day 2010. We're doing this just for fun, of course. The method used here is terribly inefficient. No doubt anyone reading this knows that google can give you a slice of pi. It's also easy to find a page with 100,000 digits of pi. Most calculators and calculator apps have a π key, of course. If you can't find a few digits of π, you're not looking very hard. There is a very nice page at Mathworld that is loaded with approximations for π. Check it out.

Ok. Let's get down to estimating the first few digits of π with a random number generator. I'll be using R and the default random number generator in R, which is the "Mersenne twister." It's used in R, python, php, ruby, and many other software systems. It has a period of $2^{19937} -1$; its "seed" is a vector of 626 integers. To be technically accurate, we should refer to it as a pseudorandom number generator (PRNG), since it is an algorithm. It is not truly random since its output is determined by initial conditions.

To calculate π, we'll focus on the first quadrant of the x,y plane within the square (0,0),  (1,0),  (1,1), and (0,1). We'll use our PRNG to generate $N = 2^{17}$ random pairs and then graph them. We'll be using a uniform random number generator, of course. Here is the R code used to generate the graph:

N=2^17
c1=6
r = function(z) { return(sqrt(z[1]^2 + z[2]^2)) }
x = runif(N) ; y = runif(N)
color = rep(c1,N)
dists = apply(cbind(x,y),1,r)
color[dists > 1] = 1
plot(y ~ x, pch=".", col = color)
ins = length(color[color == c1])
est = 4 * ins / N
cat(sprintf("estimate for pi = %.4f",est))

The apply() function is used here to calculate the distance from each point to the origin. Our interface to the PRNG is runif(), which returns a random uniform deviate. If the point lies within the unit circle, it is colored. If it lies outside, it will be black. When run, the following graph is generated:

Now, since the points were chosen with a uniform distribution, we expect a fraction π/4 of them to be colored, and 1 - π/4 to be black. The program simply counts the number of colored points and estimates π as 4 times this number divided by N. When I ran it just now, the program printed

estimate for pi = 3.1428.

So that's our Monte Carlo estimate for π.