function monte2
%-----------------------------------------------------------------------
% Normal Dist. Monte Carlo Comparison Example: Numerical Integration % 
%  (monte2 1-dim monte carlo (MCS507 Handout MCA2))      
%-----------------------------------------------------------------------
% Initialize
       clc                 % clear command window 
       clear               % clear variables
       a = -1;             % lower limit
       b = 1;              % upper limit
       n = 100;            % number of panels
       NMC = 1000000;     % select random sample size
       f = inline ('exp(-x.*x/2) ./ sqrt(2*pi)','x');      
fprintf('Monte Carlo Example: Numerical Integration Comparison\n');
       h = (b - a)/n; 
% Trapezoid Rule (see also MATLAB Trapz Function):
       tr=(f(a)+f(b))/2;
       for i = 1:n-1, tr=tr+f(a+i*h); end
       tr = h*tr;
       fprintf('(%3i-point Trapezoid: I(-1,1) = %.6f\n',n+1,tr)
% Simpson's Rule:
       sr = f(a)+f(b);
       for i = 1:n-1, if mod(i,2), sr=sr+ 4*f(a + i*h); else, sr=sr+2*f(a + i*h); end, end  
       sr = h*sr/3;
       fprintf('%3i-point Simpson: I(-1,1) = %.6f\n',n+1,sr)
% MATLAB Quad (Adaptive Simpson's Rule with Default 1.e-6 Accuracy) Function:
       qr = quad(f,a,b);
       fprintf('%8.1e-accurate Quad: = %.6f\n',1.e-6,qr)
% Direct Monte Carlo Rejection method
       fprintf('\nMonte Carlo Rejection method:\n')
       c = 0; x = randn(NMC,1);
       for k = 1:NMC
          if (x(k) >= a) & (x(k) <= b)
             c = c + 1;
          end
          if (k==10) | (k==100) | (k==1000) | (k==10000) | (k==100000) | (k==1000000)
             fprintf ('N = %6i  %.6f \n',k,c/k);
          end
       end
%----------------------------------------------------------------