# Monte Carlo method to estimate pi, but I'm lost

1 view (last 30 days)

Show older comments

Nd2 = [10, 100, 1000, 1e4, 1e5, 1e6]; % vector of sample sizes.

pihat = zeros(size(Nd2)); % vector to contain the estimates of pi at each sample size.

for i=1:length(Nd2)

???

end

% Plot figure.

figure; plot(Nd2,pihat,'bo-','LineWidth',5)

hold on; plot(Nd2,pi*ones(size(Nd2)),'r--','LineWidth',5) % true value of pi

set(gca, 'XScale', 'log') % set x axis to log scale

xlabel('$N_D$','Interpreter','latex','fontsize',20)

ylabel('$\hat{\pi}$','Interpreter','latex','fontsize',20)

title('Monte Carlo estimate of $\pi$','Interpreter','latex','fontsize',20)

So I only got as far as creating the general skeleton of my code, but the equations I put in place of ??? always comes out wrong. I'm not sure I fully understand what's going on here, so I'd appreciate it if someone can help me out. Thanks

##### 3 Comments

John D'Errico
on 4 Oct 2022

Edited: John D'Errico
on 4 Oct 2022

You've not actually done anything. :)

Your code reminds me of the old comic illustration, where a professor has some equations on a blackboard, with a cloud in the middle, that says nothing more than, "and then a miracle occurs".

Anyway, there are several schemes I can think of to estimate pi using a Monte Carlo scheme. For example, you might use a Monte Carlo scheme to estimate the area of a circle, by comparing the number of points that fall inside a unit circle, compared to the number of points that fall in a square that just contains the circle. That would give one estimate of pi, and would be trivial to write.

But another scheme might implement the old Buffon needle problem, counting the number of line segments that cross sets of parallel lines. This would be a bit more sophisticated, but it is an enirely valid scheme to estimate pi, and it is indeed a form of Monte Carlo.

Anyway, if you want more help, you should provide the actual code you tried to write. Then we could understand what variation of Monte Carlo you want to employ. Anyway, as you have it here, you are asking us to do all of the "hard" part, and that is wrong.

### Answers (1)

William Rose
on 5 Oct 2022

Edited: William Rose
on 5 Oct 2022

[correct my many typos]

The idea is to use Monte Carlo estimation to estimate pi, and to see how the estimate gets better when you use more random numbers - where Nd2 is the number of random numbers.

Nd2 = [10, 100, 1000, 1e4, 1e5, 1e6];

piEst=zeros(size(Nd2)); %allocate array for the esitmates of pi

for i=1:length(Nd2) %outer loop starts here

rand() generates a random number uniform on [0,1). So if you do

x=2*rand(Nd2(i),1)-1;

y=2*rand(Nd2(i),1)-1;

you get the x,y coordinates of Nd2 points, randomly and uniformly distributed in the -1 to +1 square. How can you use those Nd2 points to estimate pi? Hint: What fraction of the points do you expect will land inside the unit circle? Write another for loop inside the one you have already written:

numInCircle=0;

for j=1:Nd2 %inner loop starts here

if %(write a functon that is true iff x(j),y(j) is in the unit circle)

numInCircle=numInCircle+1;

end

end

Then, still inside the outer loop, you use numInCircle to estimate pi:

piEst(i)= 1; %piEst(i) should be a function of numInCircle and Nd2(i)

end %end of outer loop

Then you can plot your results.

Extra credit: Replace the inner loop with single line formula for numInCircle.

Good luck!

##### 3 Comments

William Rose
on 5 Oct 2022

Edited: William Rose
on 5 Oct 2022

[edit: correct my typos, which are not in the code]

A "then" is not needed after an "if", in Matlab. If the condition is true, the code that follows, up to the next end, will be executed. If the condition is false, execution jumps to the line after the end statement.

for i=1:5

if mod(i,2)==0

disp([num2str(i),' is even.']);

end

end

You could do that if statement on one line:

for i=1:5

if mod(i,2)==0, disp([num2str(i),' is even.']); end

end

Try it.

William Rose
on 5 Oct 2022

### See Also

### Categories

### Community Treasure Hunt

Find the treasures in MATLAB Central and discover how the community can help you!

Start Hunting!