% Monte Carlo computation of pi.
n = input(' Enter n: ');
count = 0;
% Generate random points in the square [-1,1]X[-1,1].
% The fraction of these that lie in the unit disk
% x^2+y^2 <= 1 will be approximately pi/4.
% Think of this as taking the average of N independent
% identically distributed random variables X_i, where
% X_i = 1 if point i lies in the disk, 0 otherwise.
Eofxsq = 0; % Compute expected value of X_i^2 to use in error estimate.
for i=1:n,
x = 2*rand-1; y = 2*rand-1;
if x^2 + y^2 <= 1, count = count + 1; Eofxsq = Eofxsq + 1^2; end;
end;
pi_approx = 4*(count/n),
err = pi - pi_approx,
Eofxsq = Eofxsq/n;
varx = Eofxsq - (count/n)^2; % Variance in individual approximations to pi/4.
sigx = sqrt(varx); % Std dev in individual approximations to pi/4.
sigma = 4*sigx/sqrt(n), % Std dev in total approximation to pi.
fprintf('Error should be less than %f, 68.3 percent of the time\n',sigma)
fprintf('Error should be less than %f, 95 percent of the time\n',2*sigma)