clear all
close all
% 
% Preamble
     disp('Microdamage Two Sample Problem - Permutation Test with a range of Statistics')
     lw = 3; 
     set(0, 'DefaultAxesFontSize', 16);
     fs = 15;
     msize = 10;
% 



rng(0,'twister')
% The data showing the score of normalized damage events are shown
% in the table below. There were n1 = 13 donors classified as young (<=45
% years old) and n2 = 17 classified as old (>45 years old). To calculate the
% microdamage progression score, the counts of damage events (extensions,
% surface originations, widenings, and combinations) are normalized to the
% bone area and summed up.
% Test the equality of mean scores for the two groups (young/old).
microdamage = [... 
   1     0.1582         0    0.3164    0.3164; ...
   1     0.5902    0.1180    0.2361         0; ...
   1     0.1196    0.2392    0.4785    0.1196;...
   1     0.2527    0.1264    0.6318         0;...
   1     0.2378         0    0.2378    0.2378;...
   1     0.1276         0         0    0.1276;...
   1          0         0    0.4058         0;...
   1          0         0    0.1349         0;...
   1          0         0    0.2106    0.1053;...
   1          0         0    0.1788         0;...
   1     0.5619         0    0.2809    0.4214;...
   1     0.1566         0    0.3132    0.9395;...
   1     0.6327         0    0.3163    0.2109;...
   2     0.2746         0    0.8239    0.2746;...
   2     0.3609         0    0.1203    0.1203;...
   2          0         0    0.6856    0.3428;...
   2     0.4740    0.1580    0.4740    0.1580;...
   2     0.1690         0    0.3379    0.6758;...
   2     0.4638    0.1546    0.7730    0.4638;...
   2     0.6332         0    0.7914    0.4749;...
   2     0.1622         0    0.3244         0;...
   2     0.2033         0    0.2033    0.4067;...
   2     0.4101         0    0.2734    0.1367;...
   2     0.2948         0    0.4421    0.5895;...
   2     0.6621    0.1655    0.1655    0.3311;...
   2     0.5029         0    1.0059    0.5029;...
   2     0.2053         0    0.5131    0.3079;...
   2     0.6780         0    0.3390    0.1130;...
   2          0    0.1208    0.2415    0.2415;...
   2     0.1738         0    0.2607    0.4345];
microdamage = round(1000 * microdamage)/1000;
agegroup    = microdamage(:,1);
extension   = microdamage(:,2);
surforigin  = microdamage(:,3);
widening    = microdamage(:,4);
combination = microdamage(:,5);
score = extension + surforigin +  widening + combination;

x1 = score(agegroup==1);
x2 = score(agegroup==2);
m1 = mean(x1); n1 = length(x1); s1=std(x1);
m2 = mean(x2); n2 = length(x2); s2=std(x2);

% 
sp = sqrt( ((n1-1)*s1^2 + (n2-1)*s2^2)/(n1 + n2 - 2)); %0.4460

t1 = (m1 - m2)/(sp * sqrt(1/n1 + 1/n2)); %-2.5560
t2 = (m1 - m2)/sqrt(s1^2/n1 + s2^2/n2);  %-2.5711
t3 =  m1 - m2; %-0.4200
t4 = sum(x1); % 9.5430
t5 = median(x1)-median(x2);  %-0.3400
t6 = trimmean(x1, 20) - trimmean(x2, 20); % -0.3977

df1=n1+n2-2; %28
df2 = (s1^2/n1 + s2^2/n2)^2 /((s1^2/n1)^2/(n1-1) + (s2^2/n2)^2/(n2-1)); %26.5179

pval1 = tcdf(t1, df1) %0.0082
pval2 = tcdf(t2, df2) %0.0080

% 
B=50000;
t1bs=[]; t2bs=[]; t3bs=[]; t4bs=[]; t5bs=[]; t6bs=[];

for b=1:B
allper=score(randperm(n1+n2));
x1b = allper(1:n1); x2b=allper(n1+1:n1+n2);
m1b = mean(x1b);   s1b=std(x1b);
m2b = mean(x2b);   s2b=std(x2b);
spb = sqrt( ((n1-1)*s1b^2 + (n2-1)*s2b^2)/(n1 + n2 - 2));
t1bs = [t1bs  (m1b - m2b)/(spb * sqrt(1/n1 + 1/n2))] ;
t2bs = [t2bs (m1b - m2b)/sqrt(s1b^2/n1 + s2b^2/n2)];
t3bs = [t3bs m1b - m2b];
t4bs = [t4bs sum(x1b)];
% X1bar - X2bar = Sx/n1 - Sy/n2 = Sx/n1 - (S - Sx)/n2 =
% Sx (1/n1 + 1/n2) - S/n2, S is fixed, Sx is a valid stat
t5bs = [t5bs median(x1b) - median(x2b)];
t6bs = [t6bs trimmean(x1b, 20) - trimmean(x2b, 20)]; 
end

nbin=50; hightbar = 3*B/nbin;
figure(1)
hist(t1bs, nbin)
hold on
plot([t1 t1],[0 hightbar],'r-','LineWidth', lw)
pval3 = sum( t1bs <= t1)/B      %0.0083
figure(2)
hist(t2bs, nbin)
hold on
plot([t2 t2],[0 hightbar],'r-','LineWidth', lw)
pval4 = sum( t2bs <= t2)/B      %0.0089
figure(3)
hist(t3bs, nbin)
hold on
plot([t3 t3],[0 hightbar],'r-','LineWidth', lw)
pval5 = sum( t3bs <= t3)/B      %0.0083

figure(4)
hist(t4bs, nbin)
hold on
plot([t4 t4],[0 hightbar],'r-','LineWidth', lw)
pval6 = sum( t4bs <= t4)/B      %0.0083

figure(5)
hist(t5bs, nbin)
hold on
plot([t5 t5],[0 hightbar],'r-','LineWidth', lw)
pval7 = sum( t5bs <= t5)/B      %0.0396

figure(6)
hist(t6bs, nbin)
hold on
plot([t6 t6],[0 hightbar],'r-','LineWidth', lw)
pval8 = sum( t6bs <= t6)/B      %0.0129
