%==========================================================================
% Gauss_CL_Simpson_Int.m
%==========================================================================
%
% Computes double-sided (DS) and single-sided (SS) Confidence Intervals for
% Gaussian/Normal dist'n, using simple-minded Simpson integration method.
%
% Created:          09/26/2012     13:30 hr     Prof. Steven Errede
%
% Last updated:     09/27/2012     14:35 hr {SME}
%
%==========================================================================
% Please see/read UIUC Physics 598AEM Lecture Notes 8 p. 3-8 for details:
% http://courses.physics.illinois.edu/phys598aem/598aem_lectures.html
%==========================================================================
%
%   Nsig = # sigma = (Xhi-Xo)/Xsig = (Xo-Xlo)/Xsig (Here: Xo = 0, Xsig = 1)
%   DS Prob is the area under the Gaussian within central band of +_ Nsig.
%   DS (1-Prob) is the remaining area - in the hi/lo-side tails of Gaussian.
%
%     Nsig     DS Prob(Nsig)  DS (1-Prob(Nsig))  SS Prob(Nsig)  SS (1-Prob(Nsig))
%   ========   =============  =================  =============  =================
%      0         0.000000        1.000000           0.500000       0.500000
%      1         0.682689        0.317311           0.841345       0.158655
%      2         0.954500        0.045500           0.977250       0.022750
%      3         0.997300        0.002700           0.998650       0.001350
%      4         0.999937        0.000063           0.999968       0.000032
%      5         0.999999        6x10^-7            1.000000       0.000000
%      6         1.000000        2x10^-9            1.000000       0.000000
%
%   0.674500     0.50            0.50               0.750000       0.250000
%   1.644854     0.90            0.10               0.950000       0.050000
%   1.959964     0.95            0.05               0.975000       0.025000
%   2.579964     0.99            0.01               0.995000       0.005000
%   3.290567     0.999           1x10^-3            0.999500       0.000500
%   3.890592     0.9999          1x10^-4            0.999950       0.000050
%   4.417173     0.99999         1x10^-5            0.999995       0.000005
%
%==========================================================================
clear all;
close all;

fprintf('\n');
fprintf('\n=================================================================');
fprintf('\n=== Gaussian Double-Sided & Single-Sided Confidence Intervals ===');
fprintf('\n===                Simpson Integration Method                 ===');
fprintf('\n=================================================================');

Xo   =  0.0; % True mean  of Gaussian/Normal dist'n.
Xsig =  1.0; % True sigma of Gaussian/Normal dist'n.

fprintf('\n');
fprintf('\n True mean,  Xo   = %f',Xo);
fprintf('\n True sigma, Xsig = %f',Xsig);

%======================================================================
% Change Nsig as per above table in order to obtain Prob and 1 - Prob:
%======================================================================
Nsig =  6.0; % e.g. 1.0, 2.0, 3.0, 4.0, 5.0, 6.0

fprintf('\n');
fprintf('\n # Sigma = %f',Nsig);

Npts = 2000000*Nsig*Xsig;

Xlo  =    Xo - Nsig*Xsig;
Xhi  =    Xo + Nsig*Xsig;

 dX = (Xhi - Xlo)/Npts;

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

Prob_ds  = 0.0;
X        = Xlo + 0.5*dX;
k        =   0;
kpts     =   0;
for j = 1:(Npts-1);
  F_of_X   = (1.0/(sqrt(2.0*pi)*Xsig))*exp(-(((X - Xo)^2)/(2.0*(Xsig^2))));
    
 dProb     = F_of_X*dX;
  Prob_ds  = Prob_ds + dProb;
  
  if (k == 0);
       kpts     = kpts + 1;
    %fprintf('\n j, kpts = %i %i',j,kpts);
     
       Xg(kpts) =        X;
      PDF(kpts) =   F_of_X;
      CDF(kpts) =  Prob_ds; % n.b. this is *not* the whole answer yet!
  end
  
  k        = k + 1;
  if (k > 9999); % we save only every other 10,000th set of data-points
      k =  0;
  end
    
  X        = X + dX;
end

% Now fix up the CDF(kpts):
for k = 1:kpts;
    CDF(k) = CDF(k) + 0.5*(1.0 - Prob_ds);
end

fprintf('\n');
fprintf('\n DS Prob, DS (1-Prob) = %f %f',Prob_ds,(1.0-Prob_ds));

Prob_ss = 1.0 - 0.5*(1.0-Prob_ds); % n.b. thus: Pss = 0.5*(1.0 + Pds)

fprintf('\n');
fprintf('\n SS Prob, SS (1-Prob) = %f %f',Prob_ss,(1.0-Prob_ss));

figure(01);
plot(Xg,PDF,'b');
grid on;
axis tight;
axis([-6.0 6.0 0.0 0.5]);
xlabel('x');
ylabel('f(x)');
title('Gaussian/Normal PDF f(x) vs. x');

figure(02);
plot(Xg,CDF,'b');
grid on;
axis tight;
axis([-6.0 6.0 0.0 1.0]);
xlabel('x');
ylabel('CDF(x < X)');
title('Gaussian/Normal CDF(x < X) = Int{_{-inf}^{X} f(x)*dx} vs. X');

beep();
fprintf('\n');
fprintf('\n Gauss_CL_Simpson_Int calculations finished \n');