Tuesday, November 30, 2010

Estimating Pi Using the Monte Carlo Method in MATLAB

These MATLAB functions use a large set of uniformly distributed pseudorandom ordered pairs to estimate the area of a circle and then uses the relationship A=πr2 to make an estimate of Pi.  The more random points used (i.e. the larger nsamples), the better the estimate of Pi, as shown in this plot of Pi estimate vs. sample size.  (Click to enlarge, if you're so inclined.  Notice the faint red line at y = π.)

There are two ideas I used for the MATLAB code, one using vectors and the other using only scalar variables.  The vector version seems to perform faster for large nsamples values.  More on that later if I feel like wasting more time!

Here is the code for the scalar only version:
function [MonteCarloPi] = PiEstimate(nsamples)
radius = 0.5;
radiussquared = radius^2;
inside = 0;
for i = 1:nsamples
    x = rand(1)-radius;
    y = rand(1)-radius;
    if x^2 + y^2 <= radiussquared,
        in = 1;            
    else
        in = 0;
    end
    inside = inside + in;
end
areaCircle = inside/nsamples;
MonteCarloPi = areaCircle/radiussquared

And here is the code for the vector version:
function [MonteCarloPi] = PiEstimateVector(nsamples)
radius = 0.5;
radiussquared = radius^2;
areaVector = zeros(nsamples,1);
x = rand(nsamples,1)-radius;
y = rand(nsamples,1)-radius;
rsquared = x.^2 + y.^2;
for i = 1:nsamples
    if rsquared(i) <= radiussquared,
        areaVector(i) = 1;
    end
end
areaCircle = sum(areaVector)/nsamples;
MonteCarloPi = areaCircle/radiussquared

And, finally, a video of the random points filling into the plane on the x and y intervals [-1,-1] with an inscribed circle of unit diameter.  The blue points fall inside the circle, while the red points fall outside the circle.  (Better take a seat before you watch this video.)