% ### EXintegrateMC1.m ###     10.08.14 (C. Bergevin)
% This code demonstrates integration via a two simple Monte carlo methods 
% for various functions over a specified interval  
% ** Method 1 ** - Computes Np (uniformly distributed) (x,y) pairs over the relevant range,
% than ascertains what fraction would contribute to the integral (i.e., which
% ones are 'under' the curve). 
% ** Method 2 ** - From the uniformly-distributed x-values, the function is
% evaluated and a mean such is subsequently determined (which then allows
% the definite integral to be determined)
% --> Both are compared to a trapz.m estimate (i.e., Riemann sum)
clear
% -----------------------------
% User Inputs
N= 100;      % number of 'data' points
Np= 1000;     % number of random samples to compute
intLimts= [-3 8];   % limits to evaluate integral over
% ===
% user specifies function to integrate
% generic model parameters
AA= 2;      % model parameters
BB= 2;
CC= 0.5;
% function
%z= @(x) AA./(BB-exp(-CC*x.^2));    % Gaussian
z= @(x) AA*sin(x)-3;                 % sinusoid
%z= @(x) AA*exp(-x)-BB;              % decaying exponential     
% -------------------------------------------------
x=linspace(intLimts(1),intLimts(2),N); % create x-values
y= z(x);        % determine associated y-values
% ***
figure(1); clf;     % plot the function
plot(x,y,'b-','LineWidth',2);    hold on; grid on; xlabel('x');  ylabel('y');
% ***
% Integrate via built-in Matlab function trapz.m
a1 = trapz(x,y);
msg = ['Integral calculated by trapz.m (i.e., Riemann sums)= ' num2str(a1,10)]; disp(msg);
% ***************
% Integrate via a brute-force Monte carlo routine
% 1. Create a bounding box about the function
if (min(y)>=0), yL=0;   else yL= min(y); end    % handle case if entire curve sits above x-axis
if (max(y)<=0), yH=0;   else yH= max(y); end    % handle case if entire curve sits below x-axis 
yD= yH- yL;
xL= min(x); xH= max(x); xD= xH- xL;
% 2. Create uniform distribution of points inside that box; then for each
% point determine whether it falls between the curve and x-axis (if so,
% counter is updated +1)
cntP= 0; cntN= 0;
for nn=1:Np
    xT= xD*rand(1)- abs(xL);
    yT= yD*rand(1)- abs(yL);
    zT= z(xT);  % evaluate the function at all the random x-values
    storeP(nn,:)= [xT yT];   % store away points (for visualization) 
    % consider more genral case for instances above or below x-axis
    if (yT>=0)
        if (yT<=zT), cntP= cntP+1; end
    end
    if (yT<0)  
        if (yT>=zT), cntN= cntN+1;   end
    end
    zComp(nn)= zT;  % compile all values of evaluated function
end
% estimate area as ratio determined via flagged counts
areaMC= ((cntP-cntN)/Np)*yD*xD; % diff. here accounts for area being above or below x-axis
msg = ['Integral calculated via Monte carlo Method 1 (area ratios) = ' num2str(areaMC)]; disp(msg);
plot(storeP(:,1),storeP(:,2),'r+'); title('Method 1 (area ratios)') % visualize the points used
% also report back estimate via average value
msg = ['Integral calculated via Monte carlo Method 2 (average value) = ' num2str(mean(zComp)*diff(intLimts))]; disp(msg);