% eda12_07
% generates simulated data for particle size example

dobsA=load('testA.txt');
dtrueA = 100.0; % known size of calibrations standard

% calculate useful statistics

disp('file testA.txt');

NA=length(dobsA);
disp(sprintf('N %d', NA ));

dbartrueA = dtrueA;
disp(sprintf('true mean %f', dbartrueA ));

dbarestA  = sum(dobsA)/NA;
disp(sprintf('estimated mean %f', dbarestA ));

sd2trueA=1;
disp(sprintf('true variance %f', sd2trueA ));

sdtrueA=sqrt(sd2trueA);
disp(sprintf('true root-variance %f', sdtrueA ));

sd2estA=((dobsA-dtrueA)'*(dobsA-dtrueA))/NA;
disp(sprintf('estimated variance (given true data) %f', sd2estA ));

sdestA=sqrt(sd2estA);
disp(sprintf('estimated root-variance (given true data) %f', sdestA ));

sd2estpA=((dobsA-dbarestA)'*(dobsA-dbarestA))/(NA-1);
disp(sprintf('estimated variance (given estimated mean) %f', sd2estpA ));

sdestpA=sqrt(sd2estpA);
disp(sprintf('estimated root-variance (given estimated mean) %f', sdestpA ));

ZA = abs(dbarestA-dbartrueA)/(sdtrueA/sqrt(NA));
disp(sprintf('Z %f', ZA ));

PZA = 1 - (normcdf(ZA,0,1)-normcdf(-ZA,0,1));
disp(sprintf('Probability of |Z| exceeding %f is %f', ZA, PZA ));

chi2A = ((dobsA-dtrueA)'/sdtrueA) * ((dobsA-dtrueA)/sdtrueA);
disp(sprintf('chi-squared %f', chi2A ));

Pchi2A = 1 - chi2cdf(chi2A,NA);
disp(sprintf('Probability of chi2 exceeding %f is %f', chi2A, Pchi2A ));

tA = abs(dbarestA-dbartrueA)/(sdestA/sqrt(NA));
disp(sprintf('t %f', tA ));

PtA = 1 - (tcdf(tA,NA)-tcdf(-tA,NA));
disp(sprintf('Probability of |t| exceeding %f is %f', tA, PtA ));

disp(' ');

% load second dataset
dobsB=load('testB.txt');
dtrueB=100; % known size of calibration standard

% calculate useful statistics

disp('file testB.txt');

NB=length(dobsB);
disp(sprintf('N %d', NB ));

dbartrueB = dtrueB;
disp(sprintf('true mean %f', dbartrueB ));

dbarestB  = sum(dobsB)/NB;
disp(sprintf('estimated mean %f', dbarestB ));

sd2trueB=1;
disp(sprintf('true variance %f', sd2trueB ));

sdtrueB=sqrt(sd2trueB);
disp(sprintf('true root-variance %f', sdtrueB ));

sd2estB=((dobsB-dtrueB)'*(dobsB-dtrueB))/NB;
disp(sprintf('estimated variance (given true data) %f', sd2estB ));

sdestB=sqrt(sd2estB);
disp(sprintf('estimated root-variance (given true data) %f', sdestB ));

sd2estpB=((dobsB-dbarestB)'*(dobsB-dbarestB))/(NB-1);
disp(sprintf('estimated variance (given estimated mean) %f', sd2estpB ));

sdestpB=sqrt(sd2estpB);
disp(sprintf('estimated root-variance (given estimated mean) %f', sdestpB ));

ZB = abs(dbarestB-dbartrueB)/(sdtrueB/sqrt(NB));
disp(sprintf('Z %f', ZB ));

PZB = 1 - (normcdf(ZB,0,1)-normcdf(-ZB,0,1));
disp(sprintf('Probability of |Z| exceeding %f is %f', ZB, PZB ));

chi2B = ((dobsB-dtrueB)'/sdtrueB) * ((dobsB-dtrueB)/sdtrueB);
disp(sprintf('chi-squared %f', chi2B ));

Pchi2B = 1 - chi2cdf(chi2B,NB);
disp(sprintf('Probability of chi2 exceeding %f is %f', chi2B, Pchi2B ));

tB = abs(dbarestB-dbartrueB)/(sdestB/sqrt(NB));
disp(sprintf('t %f', tB ));

PtB = 1 - (tcdf(tB,NB)-tcdf(-tB,NB));
disp(sprintf('Probability of |t| exceeding %f is %f', tB, PtB ));

disp(' ');

% are the means of the two tests significantly different ?

disp('difference between two means with Z-test');
A = dbarestA - dbarestB;
B = sqrt(sd2trueA/NA+sd2estpB/NB);
Z = abs(A/B);
disp(sprintf('Z %f', Z ));

% are the means of the two tests significantly different ?
PZ = 1 - (normcdf(Z,0,1)-normcdf(-Z,0,1));
disp(sprintf('Probability of |Z| exceeding %f is %f', Z, PZ ));
disp(' ');

disp('difference between two means with t-test');
A = dbarestA - dbarestB;
B = sqrt(sd2estpA/NA+sd2estpB/NB);
t = abs(A/B);
C = (sd2estpA/NA+sd2estpB/NB)^2;
D = ((sd2estpA/NA)^2)/(NA-1) + ((sd2estpB/NB)^2)/(NB-1);
M=round(C/D);
disp(sprintf('t %f degrees of freedom %d', t, M ));

% are the means of the two tests significantly different ?
Pt = 1 - (tcdf(t,M)-tcdf(-t,M));
disp(sprintf('Probability of |t| exceeding %f is %f', t, Pt ));
disp(' ');

% are the variances of the two tests significantly different ?
F = (chi2A/NA) / (chi2B/NB);
if( F<1 )
    F=1/F;
end
disp(sprintf('1/F %f F %f', 1/F, F ));

PF = 1 - (fcdf(F,NA,NB)-fcdf(1/F,NA,NB));
disp(sprintf('Probability of F outside interval %f %f is %f', 1/F, F, PF ));




