%=========================================================================
% MC_Exponential_Dist.m
% Written by Prof. S. Errede          last updated: 10/16/2010 10:15 hr
%
% Simple MatLab Monte Carlo exercise with the EXPONENTIAL distribution:
% a.) Generates Nevts from random Exponential distribution 
%     - only need to specify true lifetime, Ttru.
% b.) Compute (unbiased) Sample Mean,Variance & Sigma
% c.) Histograms & plots the histogram result
% d.) Converts the histogram into a (binned) P.D.F., plots the result
% e.) Calculates the variance & covariance assoc. w/ histo bin contents
%
% n.b. each time this script is run, will get different MC results - 
%      the (internal) random # seed is based on time-of-day...
%=========================================================================
close all; % clear out all figures, etc. from immediately previous use..
clear all; % clear/zero-out/initialize all variables & arrays...

int32       Nhist(100);
int32       Mhist(100);
int32       Lhist(100);

double    LnNhist(100);

double      Phist(100);
double      Chist(100);

double       Xctr(100);

double    Xrand(10000);
double    Crand(10000);
double    Prand(10000);

double   Hist_var(100);
double   Hist_cov(100,100);
double ALnHst_cov(100,100);

double Ttru;
double Vtru;
double Stru;

double Xlo;
double Xhi;

double  Xsum;
double X2sum;
double  Xavg;
double X2avg;
double  Xvar;
double  Xsig;

double Phist_sum;

% Define the # events per "experiment":
Nevts = 10000;

fprintf(' \n');
fprintf('# MC Events = %i \n',Nevts);

% MC generate Nevnts from Exponential distribution 
% with true mean Ttru .
% Type "help random <enter>" in MatLab command window for more details.
% Mouse-click on blue-highlighted "doc random" for more info...

fprintf(' \n');
fprintf('MC Exponential Distribution \n');

Ttru =  10.0;
Vtru = Ttru^2;
Stru = sqrt(Vtru);

fprintf(' \n');
fprintf('True Mean     = %f \n',Ttru);
fprintf('True Variance = %f \n',Vtru);
fprintf('True Sigma    = %f \n',Stru);

for i=1:Nevts;
    Xrand(i) = random('exp',Ttru);
end;

% Calculate the (unbiased) Sample Mean, Sample Variance and Sample Sigma:
Xsum  = 0.0;
X2sum = 0.0;
for i = 1:Nevts;
    Xsum  =  Xsum +  Xrand(i);
    X2sum = X2sum + (Xrand(i))^2;
end;

Xavg  =  Xsum/Nevts; % *unbiased* sample mean = simple/arithmetic mean
X2avg = X2sum/Nevts;
Xvar  = (Nevts/(Nevts-1))*(X2avg-(Xavg)^2); % *unbiased* sample variance
Xsig  =  sqrt(Xvar); % *unbiased* sample sigma

fprintf(' \n');
fprintf('<x>    = %f \n',Xavg); % should be ~  50.0
fprintf('Var(x) = %f \n',Xvar); % should be ~ 100.0
fprintf('Sig-x  = %f \n',Xsig); % should be ~  10.0

% Set up a 100 bin histogram on the interval [0,100]
% Define Xlo, Xhi and the bin width, bin centers:
Xlo   =   0.0;
Xhi   = 100.0;

Nbins = 100;
Xint  = Xhi - Xlo;
dX    = Xint/Nbins;

%fprintf(' \n');
for i = 1:Nbins;
    Xbin_lo = Xint*((i-1)/Nbins);
    Xctr(i) = Xbin_lo + (dX/2.0);
    
    % type/print out the Xctr array:
    % fprintf('%i %f \n',i,Xctr(i));
end;

% Now histogram the above MC generated data:
% Type "help hist <enter>" in the MatLab command window for more details:
Nhist = hist(Xrand,Xctr);

% Type/print out the Nhist array:
%fprintf(' \n');
%for i = 1:Nbins;
%    fprintf('%i %i \n',i,Nhist(i));
%end;
    
% Plot the histogrammed/binned data as a (blue) bar graph:
% type "help figure <enter>" and/or "help bar <enter>"
% in the MatLab command window for more details.
figure(01);
%plot(Xctr,Nhist,'b');
bar(Xctr,Nhist,'b');
grid on;
xlabel('x');
ylabel('N(x)');
title('Exponential Distribution: N(x) vs. x')

% Take (natural) log_e of # events in each histogram bin:
for i = 1:Nbins;
    LnNhist(i) = 0.0;
    if (Nhist(i) > 0)% can't take log(0)!
        LnNhist(i) = log(Nhist(i));
    end
end;

% Plot the log_e of histogrammed/binned data as a (blue) bar graph:
% type "help figure <enter>" and/or "help bar <enter>"
% in the MatLab command window for more details.
figure(02);
%plot(Xctr,LnNhist,'b');
bar(Xctr,LnNhist,'b');
grid on;
xlabel('x');
ylabel('Ln N(x)');
title('Exponential Distribution: Ln N(x) vs. x')

% Now we normalize the MC Nhist data to turn it into a MC PDF:
for i = 1:Nbins;
    Phist(i) = Nhist(i)/Nevts;
end;

% Type/print out the Phist array:
%fprintf(' \n');
%for i = 1:Nbins;
%    fprintf('%i %f \n',i,Phist(i));
%end
    
% Plot the PDF data with red diamonds and w/ blue joining lines:
% type "help figure <enter>" and/or "help plot <enter>" 
% in the MatLab command window for more details.
figure(03);
plot(Xctr,Phist,'rd',Xctr,Phist,'b-');
grid on;
xlabel('x');
ylabel('P(x)');
title('Exponential Distribution: P(x) vs. x')

figure(04);
semilogy(Xctr,Phist,'rd',Xctr,Phist,'b-');
grid on;
xlabel('x');
ylabel('Ln P(x) ');
title('Exponential Distribution: Ln P(x) vs. x')

% Calculate Cumulative PDF from Phist data:
Phist_sum = 0.0;
for i = 1:Nbins;
    Phist_sum = Phist_sum + Phist(i);
    Chist(i)  = Phist_sum;
end;

% Plot the Cumulative PDF data with magenta stars and w/ blue joining lines:
% type "help figure <enter>" and/or "help plot <enter>" 
% in the MatLab command window for more details.
figure(05);
plot(Xctr,Chist,'m*',Xctr,Chist,'b-');
grid on;
xlabel('x');
ylabel('C(x)');
title('Exponential Distribution: C(x) vs. x')

% Calculate the Variance and Covariance associated w/ histogram bins:
% due to multinomial probability distribution assoc/ w/ N histogram bins
for i = 1:Nbins;
    Hist_var(i)   = Nevts*Phist(i)*Phist(i);
    for j = 1:Nbins;
          Hist_cov(i,j) = -Nevts*Phist(i)*Phist(j); % n.b. (i,j) symmetric!
        AlnHst_cov(i,j) = 0.0;
        if (Hist_cov(i,j) < 0.0) % can't take log(0)!
            ALnHst_cov(i,j) =  log(abs(Hist_cov(i,j)));
        end
    end;
      Hist_cov(i,i) = 0.0; % diagonal elements of cov(x,x) don't exist
    ALnHst_cov(i,i) = 0.0; % diagonal elements of cov(x,x) don't exist
end;

% Plot the histogram variance data with green circles and w/ blue joining lines:
% type "help figure <enter>" and/or "help plot <enter>" 
% in the MatLab command window for more details.
figure(06);
plot(Xctr,Hist_var,'go',Xctr,Hist_var,'b-');
grid on;
xlabel('x');
ylabel('Hist var(x)');
title('Exponential Distribution: Histogram Bin Variance')

figure(07);
semilogy(Xctr,Hist_var,'go',Xctr,Hist_var,'b-');
grid on;
xlabel('x');
ylabel('Hist Ln var(x)');
title('Exponential Distribution: Ln of Histogram Bin Variance')

% Plot the histogram covariance data as a 3-D surface:
% type "help figure <enter>" and/or "help surf <enter>" 
% in the MatLab command window for more details.
figure(08);
surf(Xctr,Xctr,Hist_cov);
shading interp;
xlabel('x');
ylabel('x');
zlabel('Hist cov(x,x)');
title ('Exponential Distribution: Histogram Bin Covariance(x,x)');

figure(09);
surf(Xctr,Xctr,ALnHst_cov);
shading interp;
xlabel('x');
ylabel('x');
zlabel('Log[abs{Hist cov(x,x)}]');
title ('Exponential Distribution: Log[Abs{Histogram Bin Covariance(x,x)}]');

%===================================================================
% Get the Cumulative P.D.F. associated with *this* Exponential P.D.F.
for i=1:Nevts;
    Crand(i) = cdf('exp',Xrand(i),Ttru);
end;

Xlo = 0.0;
Xhi = 1.0;

fprintf(' \n');
fprintf('Xlo = %f \n',Xlo);
fprintf('Xhi = %f \n',Xhi);

% Set up a 100 bin histogram on the interval [0,1]
% Define the bin width, bin centers:
Nbins = 100;
Xint  = Xhi - Xlo;
dX    = Xint/Nbins;

%fprintf(' \n');
for i = 1:Nbins;
    Xbin_lo = Xint*((i-1)/Nbins);
    Xctr(i) = Xbin_lo + (dX/2.0);
    
    % type/print out the Xctr array:
    % fprintf('%i %f \n',i,Xctr(i));
end;

% Now histogram the above MC generated Cumulative P.D.F. data:
% Type "help hist <enter>" in the MatLab command window for more details:
Mhist = hist(Crand,Xctr);

% Type/print out the Mhist array:
%fprintf(' \n');
%for i = 1:Nbins;
%    fprintf('%i %i \n',i,Mhist(i));
%end
    
% Plot the histogrammed/binned data as a (blue) bar graph:
% type "help figure <enter>" and/or "help bar <enter>"
% in the MatLab command window for more details.
figure(10);
%plot(Xctr,Mhist,'b');
bar(Xctr,Mhist,'b');
grid on;
xlabel('x');
ylabel('N(x)');
title('Cumulative Exponential Distribution: N(x) vs. x')

% Now we normalize the MC Mhist data to turn it into a MC CDF:
for i = 1:Nbins;
    Phist(i) = Mhist(i)/Nevts;
end;

% Type/print out the Phist array:
%fprintf(' \n');
%for i = 1:Nbins;
%    fprintf('%i %f \n',i,Phist(i));
%end
    
% Plot the CDF data with red diamonds and w/ blue joining lines:
% type "help figure <enter>" and/or "help plot <enter>" 
% in the MatLab command window for more details.
figure(11);
plot(Xctr,Phist,'rd',Xctr,Phist,'b-');
grid on;
xlabel('x');
ylabel('F(x)');
title('Cumulative Exponential Distribution: F(x) vs. x')

% Now calculate the P-Value = Single-Sided Upper Confidence Level (C.L.)
% associated with *this* Exponential P.D.F.
for i=1:Nevts;
    Prand(i) = 1.0 - Crand(i);
end;

Xlo = 0.0;
Xhi = 1.0;

fprintf(' \n');
fprintf('Xlo = %f \n',Xlo);
fprintf('Xhi = %f \n',Xhi);

% Set up a 100 bin histogram on the interval [0,1]
% Define the bin width, bin centers:
Nbins = 100;
Xint  = Xhi - Xlo;
dX    = Xint/Nbins;

%fprintf(' \n');
for i = 1:Nbins;
    Xbin_lo = Xint*((i-1)/Nbins);
    Xctr(i) = Xbin_lo + (dX/2.0);
    
    % type/print out the Xctr array:
    % fprintf('%i %f \n',i,Xctr(i));
end;

% Now histogram the above MC generated P-Value = SS upper C.L. data:
% Type "help hist <enter>" in the MatLab command window for more details:
Lhist = hist(Prand,Xctr);

% Type/print out the Lhist array:
%fprintf(' \n');
%for i = 1:Nbins;
%    fprintf('%i %i \n',i,Lhist(i));
%end
    
% Plot the histogrammed/binned data as a (blue) bar graph:
% type "help figure <enter>" and/or "help bar <enter>"
% in the MatLab command window for more details.
figure(12);
%plot(Xctr,Lhist,'b');
bar(Xctr,Lhist,'b');
grid on;
xlabel('x');
ylabel('P-Value N(x)');
title('P-Value = SS Upper C.L. for Exponential Distribution: P-Value N(x) vs. x')

% Now we normalize the MC Lhist data:
for i = 1:Nbins;
    Phist(i) = Lhist(i)/Nevts;
end;

% Type/print out the Phist array:
%fprintf(' \n');
%for i = 1:Nbins;
%    fprintf('%i %f \n',i,Phist(i));
%end
    
% Plot the P-value data with red diamonds and w/ blue joining lines:
% type "help figure <enter>" and/or "help plot <enter>" 
% in the MatLab command window for more details.
figure(13);
plot(Xctr,Phist,'rd',Xctr,Phist,'b-');
grid on;
xlabel('x');
ylabel('P-value(x)');
title('P-Value = SS Upper C.L. for Exponential Distribution: P-Value(x) vs. x')

%==========================================================================
fprintf('\n MC_Exponential_Dist completed !!! \n')
%==========================================================================