(This article was originally published at The DO Loop, and syndicated at StatsBlogs.)

*#PiDay computation teaches lesson for every day: Choose estimator with smallest variance #StatWisdom*

Click To Tweet

### Monte Carlo estimates of pi

To compute Monte Carlo estimates of pi, you can use the function*f*(

*x*) = sqrt(1 –

*x*

^{2}). The graph of the function on the interval [0,1] is shown in the plot. The graph of the function forms a quarter circle of unit radius. The exact area under the curve is π / 4. There are dozens of ways to use Monte Carlo simulation to estimate pi. Two common Monte Carlo techniques are described in an easy-to-read article by David Neal (

*The College Mathematics Journal*, 1993). The first is the "average value method," which uses random points in an interval to estimate the average value of a continuous function on the interval. The second is the "area method," which enables you to estimate areas by generating a uniform sample of points and counting how many fall into a planar region.

### The average value method

In calculus you learn that the average value of a continuous function*f*on the interval [

*a, b*] is given by the following integral: In particular, for

*f*(

*x*) = sqrt(1 –

*x*

^{2}), the average value is π/4 because the integral is the area under the curve. In symbols, If you can estimate the left hand side of the equation, you can multiply the estimate by 4 to estimate pi. Recall that if

*X*is a uniformly distributed random variable on [0,1], then

*Y*=

*f*(

*X*) is a random variable on [0,1] whose mean is f

_{avg}. It is easy to estimate the mean of a random variable: you draw a random sample and compute the sample mean. The following SAS/IML program generates

*N*=10,000 uniform variates in [0,1] and uses those values to estimate f

_{avg}= E(f(X)). Multiplying that estimate by 4 gives an estimate for pi.

proc iml; call randseed(3141592); /* use digits of pi as a seed! */ N = 10000; u = randfun(N, "Uniform"); /* U ~ U(0, 1) */ Y = sqrt(1 - u##2); piEst1 = 4*mean(Y); /* average value of a function */ print piEst1; |

### The area method

Consider the same quarter circle as before. If you generate a 2-D point (an ordered pair) uniformly at random within the unit square, then the probability that the point is inside the quarter circle is equal to the ratio of the area of the quarter circle divided by the area of the unit square. That is,*P*(point inside circle) = Area(quarter circle) / Area(unit square) = π/4. It is easy to use a Monte Carlo simulation to estimate the probability

*P*: generate

*N*random points inside the unit square and count the proportion that fall in the quarter circle. The following statements continue the previous SAS/IML program:

u2 = randfun(N, "Uniform"); /* U2 ~ U(0, 1) */ isBelow = (u2 < Y); /* binary indicator variable */ piEst2 = 4*mean(isBelow); /* proportion of (u, u2) under curve */ print piEst2; |

### The variance of the Monte Carlo estimators

I confess: I experimented with many random number seeds before I found one that generated a sample for which the average value method produced a worse estimate than the area method. The truth is, the average values of the function usually give better estimates. To demonstrate this fact, the following statements generate 1,000 estimates for each method. For each set of estimates, the mean (called the Monte Carlo estimate) and the standard deviation (called the Monte Carlo standard error) are computed and displayed:/* Use many Monte Carlo simulations to estimate the variance of each method */ NumSamples = 1000; pi = j(NumSamples, 2); do i = 1 to NumSamples; call randgen(u, "Uniform"); /* U ~ U(0, 1) */ call randgen(u2, "Uniform"); /* U2 ~ U(0, 1) */ Y = sqrt(1 - u##2); pi[i,1] = 4*mean(Y); /* Method 1: Avg function value */ pi[i,2] = 4*mean(u2 < Y); /* Method 2: proportion under curve */ end; MCEst = mean(pi); /* mean of estimates = MC est */ StdDev = std(pi); /* std dev of estimates = MC std error */ print (MCEst//StdDev)[label="Comparison of Methods" colname={"Avg Func" "Area"} rowname={"MC Estimate" "MC StdErr"}]; |

tags: pi day, Simulation, Statistical Thinking

The post Monte Carlo estimates of pi and an important statistical lesson appeared first on The DO Loop.

**Please comment on the article here:** **The DO Loop**