%=========================================================================
% Quotient_Error_Propagation.m
% Written by Prof. S. Errede 09/17/2012 12:55 hr
%
% Last Updated:              09/17/2012 14:45 hr {SME}
%=========================================================================
close all; % clear out all figures, etc. from immediately previous use..
clear all; % clear/zero-out/initialize all variables & arrays...

Npts = 7;

Sig0 = zeros(1,Npts);
Tsig = zeros(1,Npts);

Sig0(1)  = 0.01;
Sig0(2)  = 0.02;
Sig0(3)  = 0.05;
Sig0(4)  = 0.10;
Sig0(5)  = 0.20;
Sig0(6)  = 0.50;
Sig0(7)  = 1.00;
%Sig0(8)  = 2.00;
%Sig0(9)  = 5.00;
%Sig0(10) = 10.0;

Xavg  = zeros(1,Npts);
Xavg2 = zeros(1,Npts);
Xvar  = zeros(1,Npts);
Xsig  = zeros(1,Npts);

Yavg  = zeros(1,Npts);
Yavg2 = zeros(1,Npts);
Yvar  = zeros(1,Npts);
Ysig  = zeros(1,Npts);

Zavg  = zeros(1,Npts);
Zavg2 = zeros(1,Npts);
Zvar  = zeros(1,Npts);
Zsig  = zeros(1,Npts);

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

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

% MC generate Nevnts of random X, Y from Gaussian distribution 
% each with true mean Xtru, Ytru and sigma SigX, SigY.
% Type "help random <enter>" in MatLab command window for info/details...
% Mouse-click on blue-highlighted "doc random" for additional info/details...

fprintf('\n');
fprintf('\n Error Propagation for Z = X/Y Quotient');

for k = 1:Npts;
    fprintf('\n');
    fprintf('\n============================');
    fprintf('\n k, Sig0 = %i %f',k,Sig0(k));
    fprintf('\n============================');
    
    Xtru =  1.0;    % true mean of X
    SigX = Sig0(k); % true sigma of X
    VarX = SigX^2;  % true variance of X

    Ytru =  1.0;    % true mean of Y
    SigY = Sig0(k); % true sigma of Y
    VarY = SigY^2;  % true variance of Y

    fprintf(' \n');
    fprintf('True X-Mean     = %f \n',Xtru);
    fprintf('True X-Variance = %f \n',VarX);
    fprintf('True X-Sigma    = %f \n',SigX);

    fprintf(' \n');
    fprintf('True Y-Mean     = %f \n',Ytru);
    fprintf('True Y-Variance = %f \n',VarY);
    fprintf('True Y-Sigma    = %f \n',SigY);

    Xrnd = zeros(1,Nevts);
    Yrnd = zeros(1,Nevts);
    Zrnd = zeros(1,Nevts);

    for j=1:Nevts;
        Xrnd(j) = random('norm',Xtru,SigX);
        Yrnd(j) = random('norm',Ytru,SigY);
        Zrnd(j) = Xrnd(j)/Yrnd(j);
    end;

    % Calculate the (unbiased) Sample Mean, Sample Variance and Sample Sigma:
    Xsum  = 0.0;
    Xsum2 = 0.0;
    Ysum  = 0.0;
    Ysum2 = 0.0;
    Zsum  = 0.0;
    Zsum2 = 0.0;
    for j = 1:Nevts;
        Xsum  = Xsum  +  Xrnd(j);
        Xsum2 = Xsum2 + (Xrnd(j))^2;
    
        Ysum  = Ysum  +  Yrnd(j);
        Ysum2 = Ysum2 + (Yrnd(j))^2;
    
        Zsum  = Zsum  +  Zrnd(j);
        Zsum2 = Zsum2 + (Zrnd(j))^2;
    end;

    Xavg(k)  =  Xsum /Nevts;    % *unbiased* sample mean = simple/arithmetic X-mean
    Xavg2(k) =  Xsum2/Nevts;

    Xvar(k)  = (Nevts/(Nevts-1))*(Xavg2(k)-(Xavg(k))^2); % *unbiased* sample X-variance
    Xsig(k)  =    sqrt(Xvar(k));                         % *unbiased* sample X-sigma

    Yavg(k)  =  Ysum /Nevts;    % *unbiased* sample mean = simple/arithmetic Y-mean
    Yavg2(k) =  Ysum2/Nevts;

    Yvar(k)  = (Nevts/(Nevts-1))*(Yavg2(k)-(Yavg(k))^2); % *unbiased* sample Y-variance
    Ysig(k)  =    sqrt(Yvar(k));                         % *unbiased* sample Y-sigma

    Zavg(k)  =  Zsum /Nevts;    % *unbiased* sample mean = simple/arithmetic Z-mean
    Zavg2(k) =  Zsum2/Nevts;

    Zvar(k)  = (Nevts/(Nevts-1))*(Zavg2(k)-(Zavg(k))^2); % *unbiased* sample Z-variance
    Zsig(k)  =    sqrt(Zvar(k));                         % *unbiased* sample Z-sigma

    fprintf(' \n');
    fprintf('<x>    = %f \n',Xavg(k)); % should be ~ Xtru
    fprintf('Var(x) = %f \n',Xvar(k)); % should be ~ VarX
    fprintf('Sig-x  = %f \n',Xsig(k)); % should be ~ SigX

    fprintf(' \n');
    fprintf('<y>    = %f \n',Yavg(k)); % should be ~ Ytru
    fprintf('Var(y) = %f \n',Yvar(k)); % should be ~ VarY
    fprintf('Sig-y  = %f \n',Ysig(k)); % should be ~ SigY

    fprintf(' \n');
    fprintf('<z>    = %f \n',Zavg(k)); % should be ~ Ztru
    fprintf('Var(z) = %f \n',Zvar(k)); % should be ~ VarZ
    fprintf('Sig-z  = %f \n',Zsig(k)); % should be ~ SigZ
end

% make a thy pred for absence of the non-linear term - use the first 3 
% Zsig terms as ~ reasonably accurate approximation for power law rel'n:
Tsig(1)  =    1.0*Zsig(1);
Tsig(2)  =    1.0*Zsig(2);
Tsig(3)  =    1.0*Zsig(3);
Tsig(4)  =   10.0*Zsig(1);
Tsig(5)  =   10.0*Zsig(2);
Tsig(6)  =   10.0*Zsig(3);
Tsig(7)  =  100.0*Zsig(1);
%Tsig(8)  =  100.0*Zsig(2);
%Tsig(9)  =  100.0*Zsig(3);
%Tsig(10) = 1000.0*Zsig(1);

figure(01);
loglog(Sig0,Zsig,'b',Sig0,Zsig,'m*',Sig0,Tsig,'r');
axis tight;
grid on;
xlabel('Sig0');
ylabel('Sigma Z');
title('Sigma Z vs. Sig0 for Z = X/Y, SigX = SigY = Sig0');

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