disp('Illustration of two-sample Kolmogorov-Smirnov test:');
disp(' ');
disp('Data from textbook Table 4.1 p110 (2nd edn)');
disp('Water Diffusion across human chorioaminon tissue');
disp('  X: At term (approx 39weeks, Control)    Y: 12-26 weeks');

xx = [ 0.80  0.83  1.89  1.04  1.45  1.38  1.91  1.64  0.73  1.46 ];
mm=length(xx);
xbar=sum(xx)/mm;

disp(' ');
fprintf('Sample 1 (Xs, m=%d, Xbar=%g)\n', mm,xbar);
for i=1:mm;  fprintf('  %5.2f',xx(i));  end;  fprintf('\n');

yy = [ 1.15  0.88  0.90  0.74  1.21 ];
nn=length(yy);
ybar=sum(yy)/nn;

disp(' ');
fprintf('Sample 2 (Ys, n=%d, Ybar=%g)\n', nn,ybar);
for i=1:nn;  fprintf('  %5.2f',yy(i));  end;  fprintf('\n');

disp(' ');  disp('After sorting each sample:');
disp('Sample 1');
xx=sort(xx);
for i=1:mm;  fprintf('  %5.2f',xx(i));  end;  fprintf('\n');
disp('Sample 2');
yy=sort(yy);
for i=1:nn;  fprintf('  %5.2f',yy(i));  end;  fprintf('\n');

disp(' ');
disp('Samples sorted together (3rd column is Group: X=0, Y=1)');
zz = [xx'  zeros(mm,1);  yy' ones(nn,1)];
zz = sortrows(zz,[1 2]);
zz = [ [1:mm+nn]'  zz];    %'

disp('Columns are Ordinal,Value,Group');
fprintf('   %2d    %5.2f    %d\n', zz')  %'

disp('  ');
disp('Calculation of maximum absolute difference');
disp('   of empirical distribution functions:');
disp('     J2 = max_t |F_{m,X}(t) - F_{n,Y}(t)|');

xadd = 1/mm;  yadd=-1/nn;

disp('  ');
fprintf('Increment for X: xadd=%.2f  for Y: yadd=%.2f\n', xadd,yadd);
lmm=lcm(mm,nn);
xad2=lmm*xadd;  yad2=lmm*yadd;
fprintf('Multiplied by LCM=%d:  xad2=%d  for yad2=%df\n', lmm,xad2,yad2);

disp(' ');
disp('Table used to calculate J = LCM*J2:  Columns are');
disp('  Ordinal,Value,Group,  Term,Diffsum,Maxdif, New?');
dsum=0;  maxdif=0;
% Values and group numbers
vals=zz(:,2);  group=zz(:,3);
for k=1:mm+nn;
  if group(k)==0 jval=xad2;  else jval=yad2;  end;
  dsum=dsum+jval;
  if abs(dsum) > maxdif  maxdif=abs(dsum);  newval=1;
  else  newval=0;  end;
  fprintf ('   %2d   %5.2f  %d     %2d    %2d   %2d   %d\n',k,vals(k),group(k),jval,dsum,maxdif,newval);
  end;

disp('  ');
disp('Maxdif=JJ for Table A10, J2=J/LCM = Maxdif of empirical distributions,');
disp('  Jstar for large-sample approximation:');
JJ = maxdif;                 %% For Table A10
J2 = maxdif/lmm;             %% Max diff of empirical d.f.s
Jstar = sqrt(mm*nn/(mm+nn))*J2;  %% For large-sample approximation
JJ_J2_Jstar = [JJ  J2   Jstar]

disp('P-value for large-sample approximation by summing an infinite series:');
qsum=0;  sgn=1;  ss=Jstar;
for k=1:100000;
  eterm=exp(-2*k*k*ss*ss);
  qsum = qsum + 2*sgn*eterm;  sgn=-sgn;
  if eterm<1e-8 break;  end;  end;
fprintf ('Infinite sum stopped after %d terms\n',k);
fprintf ('(Intrinsically two-sided) approximate P-value = %6.4f\n', qsum);

disp(' ');
disp('Doing nsims permutations to estimate the P-value:');
nsims=10000

rand('state',0);
Ngr=0;
for ns=1:nsims;
    kk=randperm(mm+nn);
    rdsum=0;  rmaxdif=0;
    % Calculate J statistic for this permutatation
    for j=1:mm+nn;
      if group(kk(j))==0  rdsum=rdsum+xad2;
      else rdsum=rdsum+yad2;  end;
      if abs(rdsum)>rmaxdif rmaxdif=abs(rdsum);  end;  end;
    % See if it is larger than observed D
    if rmaxdif >= maxdif  Ngr=Ngr+1;  end;  end;

disp('Number greater-than-or-equal and number of simulations:');
disp(' ');
disp([Ngr nsims]);

phat=Ngr/nsims;  pwid=1.960*sqrt(phat*(1-phat)/nsims);
disp('95% confidence interval for simulated 2-sided P-value');
disp('  bracketed by estimated P-value:');
disp(' ');
disp([ phat-pwid  phat  phat+pwid ]);

disp(' ');
disp('Compare with result of Wilcoxon Rank-Sum Test:');
disp('Data with midranks:  Columns are');
disp('  Ordinal,Value,Group, Midrank,Tiegroups');

[midr tiegr] = getmidranks(vals);
disp([zz  midr  tiegr]);

Wx=sum(midr.*(1-group));  Wy=sum(midr.*group);
ExpWx=mm*(mm+nn+1)/2;  ExpWy=mm*(mm+nn+1)/2;
Wx_Wy_Sum_ExpWx_ExpWy = [Wx  Wy  Wx+Wy  ExpWx  ExpWy]

Wvar = mm*nn*(mm+nn+1)/12;
Zval = (Wx - ExpWx)/sqrt(Wvar);

Wilcoxon_RankSum_Zval_NormPval = [Zval  normpval(Zval)]


