function mcm2acceptreject09
% Fig. 9.10 mcm2acceptreject: Monte Carlo Method (3/2007),
%     nx = 1 dim, normal dist., 
%     I = int(h(x),x=a..b), h(x) = f(x) = exp(-x^2/2)/sqrt(2pi), -2 <= a < b <= +2;
%     technically, f(x) = I_{x in [a,b]} to account for truncated integral,  
%     so I = meanf.
%
clc; clear
%
fprintf('Monte Carlo and Finite Difference Comparison:'); 
fprintf('\n   including Acceptance-Rejection Technique Application,');
fprintf('\n   with Normal Dist. on (a,b)'); 
fprintf('\n   and with int(f(x),x=a..b), f(x) = exp(-x^2/2)/sqrt(2pi);\n');
%
a = -2; b = 2; % limits of integration;
nfd = 100; % number of finite difference steps;
kmax = 7;
nmc = 10^kmax; % select Monte Carlo random sample size;
f = inline ('exp(-x.*x/2)./sqrt(2*pi)','x');  % x in [a,b]  
%   Thus, relative to the normal density, f(x)={1, x in [a,b]; 0, else};
h = (b - a)/nfd; % step size;
% Trapezoid Rule (see also MATLAB trapz(x,y) built-in function):
trap = (f(a)+f(b))/2;
for i = 1:nfd-1, 
    trap = trap+f(a+i*h); 
end
trap = h*trap;
fprintf('\n%3i-point Trapezoidal Rule: I(-1,1) = %.6f\n',nfd+1,trap);
% Simpson's (1/3) Rule:
simp = f(a)+f(b);
for i = 1:nfd-1 
    if mod(i,2) 
        simp =simp+ 4*f(a + i*h); 
    else 
        simp=simp+2*f(a + i*h); 
    end 
end  
simp = h*simp/3;
fprintf('\n%3i-point Simpson''s rule: I(-1,1) = %.6f\n',nfd+1,simp);
% MATLAB quad built-in function (adaptive Simpson's rule, default 1.e-6 accuracy):
tol = 1.e-9;
quadfn = quad(f,a,b,tol);
fprintf('\n%7.1e-accurate quad: = %.6f\n',tol,quadfn);
% Direct von Neumann Acceptance-Rejection Technique:
fprintf('\nMonte Carlo results by von Neumann''s Acceptance-Rejection technique:\n');
fprintf('\n k          n    Ihatn    stderrn\n');
nac = 0; 
x = randn(1,nmc); % MATLAB nmc X 1 normal distribution;
kv = zeros(1,kmax); Ihatn = zeros(1,kmax); stderrn = zeros(1,kmax); 
for n = 1:nmc
    if (x(n) >= a) && (x(n) <= b)
        nac = nac + 1; % counts accepted points;
    end
    if (n==10)||(n==100)||(n==1000)||(n==10000)||(n==100000)||(n==1000000)||(n==nmc)
        k = log10(n);
        kv(k) = k;
        Ihatn(k) = nac/n; % Estimate of the re-normalization constant/Integral
        stderrn(k) = sqrt(Ihatn(k)*(1-Ihatn(k))/(n-1));
        fprintf ('%2i   %8i  %8.6f  %9.3e\n',k,n,Ihatn(k),stderrn(k));
    end
end
fprintf('\n 101-pt. trap: %8.6f  %9.3e*',trap,abs(trap-quadfn));
fprintf('\n 101-pt. simp: %8.6f  %9.3e*',simp,abs(simp-quadfn));
fprintf('\n    accurate:  %8.6f  %9.3e*',quadfn,abs(quadfn-quadfn));
fprintf('\n * Absolute Errors\n');
%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
scrsize = get(0,'ScreenSize');
ss = [3.0,2.8,2.6,2.4,2.2,2.0];
%
nfig = 1;
figure(nfig);
kv = 1:kmax;
plot(kv,Ihatn,'k-o','linewidth',3,'MarkerSize',12); hold on
plot(kv,10*stderrn,'k-x','linewidth',3,'MarkerSize',12); hold off
axis([min(kv) max(kv) 0 1]);
title('Monte Carlo, Normal Dist., h(x) = f_n(x) on [a,b]'...
    ,'Fontsize',44,'FontWeight','Bold');
xlabel('log(n), Log_{10} Sample Size'...
    ,'Fontsize',44,'FontWeight','Bold');
ylabel('Moments I_n, 10*std-err_n'...
    ,'Fontsize',44,'FontWeight','Bold');
legend('I_n, Integral-est.','10*std-err_n','Location','Best');
set(gca,'Fontsize',36,'FontWeight','Bold','linewidth',3);
set(gcf,'Color','White','Position' ...
    ,[scrsize(3)/ss(nfig) 70 scrsize(3)*0.60 scrsize(4)*0.80]);
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%
% end mcm2acceptreject.m
%