Using loop to approximate pi (Monte Carlo, MATLAB)

Question

I've written the following code, based on a for loop to approximate the number pi using the Monte-Carlo-method for 100, 1000, 10000 and 100000 random points.

count = 0;      % count variable, start value set to zero
n(1)=100;
n(2)=1000;
n(3)=10000;
n(4)=100000;
for k=1:4;
  for i=1:n(k);   % for loop initialized at one
    x=rand;     % rand variable within [0,1] ... points coordinates
    y=rand;     

    if (x^2+y^2 <=1)  %if pair is within the first quadrant of the unit circle
    count = count +1;   % if yes, it increases the count by one
    end
  end
end

piapprox = 4*(count/n(k));

The whole code makes sense to me but it seems to be wrong. I couldn't finish it so far. Could you tell me please what exactly is wrong?

Answer 1

There are some issues in the code. Here is a revised version:

n(1)=100;
n(2)=1000;
n(3)=10000;
n(4)=100000;
piapprox = zeros(1, 4); % Save approximations
for k=1:4;
  count = 0;      % count variable, start value set to zero
  for i=1:n(k);     % for loop initialized at one
    x=rand;         % rand variable within [0,1] ... points coordinates
    y=rand;     

    if (x^2+y^2 <=1)  %if pair is within the first quadrant of the unit circle
    count = count +1;   % if yes, it increases the count by one
    end
  end
  piapprox(k) = count;
end

piapprox = 4 .* (piapprox ./ n);

• You need to reset the counter for every $k$
• You need to save the individual counts or the first three iteraitons of the $k$ loop will be useless

Answer 2

This should do the job:

n(1)=100;
n(2)=1000;
n(3)=10000;
n(4)=100000000;
for k=1:4;
count = 0; % reset count to zero for each test
for i=1:n(k); % for loop initialized at one
x=rand; % rand variable within [0,1] ... points coordinates
y=rand;

if (x^2+y^2 <=1)  
     count = count +1;  
end

end
piapprox(k) = 4*(count/n(k)); % keep all the values of pi approxmated
end

Answer 3

In Matlab you could easily vectorize the inner loop. That will make the program neater, and will reduce running time:

n = [100 1000 10000 100000];
piapprox = NaN(size(n)); %// initiallize result
for k = 1:numel(n)
    piapprox(k) = 4*sum(sum(rand(2,n(k)).^2) < 1)/n(k);
end
disp(piapprox)