% mc.m
% Monte Carlo simulation to compute the value of pi
% by finding the fraction of points lying inside the
% unit circle within the first quadrant.
% K. Ming Leung, September 17, 2008.

clear;

% Running only one simulation at a time:
disp('Running only one simulation at a time:');
n0 = 10000;  % total # of pts
Xvec = rand(1,n0);
Yvec = rand(1,n0);

n = sum(Xvec.^2 + Yvec.^2 <= 1);
myPi = 4 * n / n0,
relativeError = myPi/pi - 1,


% Running a several simulations at a time:
disp('==========================================');
disp('Running several simulations at a time:');
nTotal = 10000;    % total # of pts
repeat = 7;     % number of simulation copies 
XX = rand(repeat,n0);   % repeat x n0 dimensional matrix
YY = rand(repeat,n0);   % repeat x n0 dimensional matrix

NumberPtsInCircle = sum(XX.^2 + YY.^2 <= 1,2);
MyPi = 4 * NumberPtsInCircle / nTotal,
RelativeErrors = MyPi/pi - 1,
myFinalPi = mean(MyPi), % take the mean for final value
finalRelError = myFinalPi/pi - 1,
