piEstimation.m

% x represents series of number of random numbers
x = 100:2000;
% y represents estimated value of PI
y = zeros(1, length(x));

for i=1:length(x)
    inliers = 0;
    % get single value from the series
    N = x(i);
    % generate N random 2-D coordinates
    coords = rand(N, 2);
    % move the uniform [0, 1] range to uniform [-1, 1] 
    randValues = -1.0 + 2 .* coords;
    for j=1:size(randValues, 1)
        % check if the point is in unit circle and count the numbers of
        % inliers
        coord = randValues(j, :);
        dist = norm(coord);
        if dist <= 1
            inliers = inliers + 1;
        end
    end
    % probability of being inlier
    % prob = area of circle / area of square
    % prob = pi * r * r / (2r * 2r) = pi / 4
    % so, pi = prob * 4
    prob = (inliers / N );
    est_pi = prob * 4;
    fprintf('Estimated PI for N = %d -> %f\n', N, est_pi) 
    y(i) = est_pi;
end

 %set up
average_pi = sum(y) / length(y);
vec_avg = average_pi .* ones(1, length(y));
fprintf('Average PI in the long run: %f\n', average_pi);

% plotting
figure
plot(x, y, 'r')
title('PI Estimation')
xlabel('Number of random numbers')
ylabel('Estimated valıue of PI')
hold on
plot(x, vec_avg, 'b', 'LineWidth', 2)