[Math] Using loop to approximate pi (Monte Carlo, MATLAB)

estimationMATLABmonte carlopirandom

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?

Best Answer

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)