"Global pond drying"

"The first step to recovery is admitting you have a problem."

"Supportive criticism"

Monte Carlo integration: How the number of iterations used affects precision

While it may be obvious that with Monte Carlo integration, the more iterations performed will result in higher precision, it's always fun to look at something graphically. Since I had already written a program that approximated the value of pi (check it out here), I decided to look at how standard deviations change with increasing numbers of iterations.

I chose to run 10 sets of iterations: 2, 4, 8, 32, 64, 128, 256, 512, and 1024. Why these values were chosen becomes obvious in the script, but you can use any set of iterations you want. Each set of iterations was run 10,000 times and standard deviations were calculated from these.

Figure 1: Stdev vs. number of iterations. Note the x-axis is log scale so that the trend can more easily be seen.

Here's my script:

% Graphing standard deviations of pi vs. number of iterations

% using Monte Carlo integration

close all

clear all

for x = 1:10;           

   

y=2^x;

   

for i = 1:10000;

    n = 2*(rand(2,y))- 1;    % generates random numbers between -1 and 1

   

    pts = 0;                  % starting value for total points

    circ = 0;                 % starting value for points within circle

   

    for R = 1:y;

        d = sqrt((n(1,R))^2 + (n(2,R))^2);

        pts = pts + 1;

        if d <= 1;

            circ = circ + 1;

        end  

    end

    Pi(1,i) = (4*circ)/pts;

end

STD(1,x) = std(Pi);

IT(1,x) = y;

end

plot (log10(IT),STD)

title('Standard Deviation vs. Number of Iterations');

xlabel('log10(number of iterations)');

ylabel('Standard deviation');

Calculating Pi: Monte Carlo Integration

Thank you physicists, Matlab, and (pseudo)random number generators. Welcome to Monte Carlo integration. A couple of days ago, my boyfriend was telling me how he has surpassed me in his Matlab skills. Dorky, I know.  Not one to be bettered, I set out to prove him wrong. He had an assignment to approximate the value of pi using Monte Carlo integration in C. I'd do it in Matlab.

But I had to figure out what Monte Carlo integration was first.

Enter Wikipedia. I don't care what professors say about it, Wikipedia is a magical place filled with more knowledge than could ever be stuffed into my head. I'm not going to start citing the website in anything important, but I will give it my thanks. So here it is: Monte Carlo integration offers an approximate evaluation of definite integrals using random numbers.   

So what did I have to do exactly? My bf gave me a hint: use a circle and a square. After some more prodding he added, "put the circle in the square." Alright. Thanks a lot.

So I drew up a circle with radius = 1, inside a square with sides = 2 (Figure 1).

Figure 1: Circle with radius = 1 within a square with sides = 2.

We know that the area of the square is L²= 2² = 4.

We know the area of the circle is
A = πr². To find π, we need to know the area of the circle. Unfortunately, that's impossible to know without π.

But there is another way to find pi. I first set the center of the circle as (0,0) in a coordinate plane. I then could generate sets of random numbers between -1 and 1 and use these as x and y coordinates. Using Pythagorean's Theorem I can find the distance of these points from the center (the origin):

c = √(a² + b²).

If a point's distance from the origin is < 1, then it falls within the circle. The ratio of points within the circle to the total number of points (points within the square) is proportionate to the ratio of the area of the circle to the area of the square:

                                               # points in circle   =    πr²           where r=1

                                                # points total             4

So:    

                                              π    =    4(# points in circle)

                                                              # points total

With that in mind, my Matlab adventure began.

The script was pretty easy to write (find it at the bottom of this page). The hardest part for me was actually plotting the figure. Obviously, the more iterations (# of points) you do, the more accurate your result will be (Table 1). It also vastly increases the computing time. Past 1,000,000 iterations, my Matlab program got pissed and said it didn't have enough memory.

# of Iterations	Approximate value of Pi
10	2.4
100	3.24
1000	3.092
10000	3.1104
100000	3.13168
1000000	3.143080

                            

Here's my Matlab script for 1000 iterations:

% Calculating pi using Monte Carlo integration

close all

clear all

n = 2*(rand(2,1000))- 1;        % generates random numbers between -1 and 1

t=sum(n.*n,1)<1;

pts = 0;                       % starting value for total points

circ = 0;                      % starting value for points within circle

for R = 1:1000;

    d = sqrt((n(1,R))^2 + (n(2,R))^2);

    pts = pts + 1;

    if d <= 1;

        circ = circ + 1;

    end

end

format long

Pi = (4*circ)/pts

plot(n(1,~t),n(2,~t),'b.','MarkerSize',15)

hold on

plot(n(1,t),n(2,t),'r.','MarkerSize',15)

theta = linspace(0,2*pi,100);    

x = cos(theta);                 

y = sin(theta);                  

plot(x,y,'k'); 

axis ('equal');

axis([-1 1 -1 1]);

title('100 Iterations');