Loren on the Art of MATLAB

Turn ideas into MATLAB

Estimating pi Using Buffon’s Method 6

Posted by Loren Shure,

I recently attended the ICIAM meeting in Valencia, Spain which meant I got to hang out with my pals Carlos Sanchis and Lucas Garcia :-)! Carlos showed me a problem he was working with Professor Fernando Giménez from UPV regarding an app for estimating $\pi$ using Buffon's method. Here's the problem statement from Wikipedia:

Suppose we have a floor made of parallel strips of wood, each the same width, and we drop a needle onto the floor. What is the probability that the needle will lie across a line between two strips?

Interesting that the original intention had nothing to do with computing $\pi$ ! There's some fun, powerful, yet fairly easy code to demonstrate the algorithm.

Contents

Set Up Parameters

How many line segments?

N = 1000;

Length of each line?

L = 0.20;

We want the beginning points of the lines to lie between L and 1-L so we don't go outside the unit square.

xb = L + rand(1,N)*(1-2*L);
yb = L + rand(1,N)*(1-2*L);
angs = rand(1,N)*360;
xe = xb + L*cosd(angs);
ye = yb + L*sind(angs);

Visualize the Lines

ax = axes;
plot(ax,[xb;xe],[yb;ye])
axis square

Show the Vertical Grid Lines Defined by L Spacing

hold on
glines = 0:L:1;
for i = 1:length(glines)
   xline(ax, glines(i));
end

Count the Segments Intersecting the Grid

n = sum(floor(xb/L) ~= floor(xe/L));
piEstimate = 2 * N / n
piEstimate =
    3.1153

Annotate Final Plot

title("Estimate of \pi is " + piEstimate)

What Happens as L and N change?

This could be a great exercise for the classroom - seeing how the estimates depend on how many line segments and the spacing of the grid. Not to mention running a bunch of times with different random numbers each time. What simple estimation problems do you like to use? Let me know here.


Get the MATLAB code

Published with MATLAB® R2019a

6 CommentsOldest to Newest

Cris Luengo replied on : 1 of 6
You’ve introduced a bias by having all needles start on the right of the leftmost grid line, and on the left of the rightmost one. This experiment is best performed with a single grid line, producing needles that fall within reach. Though I think it’s silly to use cos and sin to get an estimate of pi... :)
@Cris- You are right - we are missing possibilities this way, e.g., going across the corners at a diagonal... I'll have to think more about that. --loren
Cris Luengo replied on : 3 of 6
Loren, I was wrong, it looks like your code doesn't produce a biased estimate. Sorry for the knee-jerk reaction there. :) This is the same concept, simplified to some extreme. It is equivalent to generating needles of length 2 by selecting a center point in the range x=[-1,1], where it could intersect the line at x=0. Note that the y-axis location is irrelevant!
x = rand(1,N) * 2 - 1;
angs = rand(1,N)*360;
n = sum(abs(x) < abs(cosd(angs)));
results(ii,2) = 2 * N / n;
Daniel M replied on : 5 of 6
I was also not very confortable using cosd and sind to calculate pi ;) but you can replace it with:
N_rand = 1000;
x_rand = rand(1, N_rand) - 0.5;
y_rand = rand(1, N_rand) - 0.5;
r_rand = sqrt(x_rand.^2 + y_rand.^2);
ix = r_rand <= 0.5;
cos_rand = x_rand(ix)./r_rand(ix);
sin_rand = y_rand(ix)./r_rand(ix);
N = sum(ix);
The only problem that you will not control N. Daniel
@Daniel- Thanks for the snippet. You can control N by generating a large enough collection by increasing N_rand - but it's not deterministic hope my tries. Just keep adding more attempts until you reach N suitable ones. --loren

Add A Comment

Your email address will not be published. Required fields are marked *

Preview: hide