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 π.

## No comments:

## Post a Comment