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

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

[Image Analyst]

What are the equations? (Don't make us look them up, because we probably won't.)

[John D'Errico]

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.

https://en.wikipedia.org/wiki/Buffon%27s_needle_problem

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.

[Kevin Nelson]

This is what I put in place of the ???, it just didn't make much sense to me though

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

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

r = 1;

k = sqrt((x-5).^2+(y+5).^2)<=r;

k(k==0) = nan;

ratio=(sum(k>0)/Nd2);

pi_value = ratio*4;

pihat(i) = pi_value;

Answer 1

[William Rose]

@Kevin Nelson,

[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!

[Kevin Nelson]

I'm confused on the if line that you wrote. Shouldn't there be a then statement that follows it. Otherwise, the code won't work

[William Rose]

@Kevin Nelson,

[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

2 is even.
4 is even.

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

2 is even.
4 is even.

Try it.

[William Rose]

@Kevin Nelson, I had not looked carefully at your reply to the question from @Image Analyst. Now that I have read your reply, I can see that the code, which did not make sense to you, was an unsuccessful attempt to estimate pi by the same approach I recommended: using two random numbers to generate a point on the x-y plane, many times. Count how many of those points are inside the unit circle. Use the result to estimate pi. I agree with you: the code does not quite make sense. But now you know what that code was trying to do.