%
% Large-sample approximation and simulation of true P-values
%   for Wilcoxon rank-sum test
% Use  getmidranks  function to find midranks

echo off;

disp('Water Diffusion across human chorioaminon tissue');
disp('  X: At term (approx 39weeks, Control)    Y: 12-26 weeks');
disp(' ');    % Blank line in output
disp('Data from textbook Table 4.1 p110 (2nd edn)');
disp('Use Wilcoxon Rank-Sum-test procedures for two samples:');

xx = [ 0.80  0.83  1.89  1.04  1.45  1.38  1.91  1.64  0.73  1.46 ];
yy = [ 1.15  0.88  0.90  0.74  1.21 ];

nx=length(xx);    % Sample sizes
ny=length(yy);
%% Size of two samples together
nn=nx+ny;

% Display the two samples
% (There is no easy way to do this in MATLAB other than by using fprintf)
fprintf('\nX (m=%2d):',nx);
for i=1:nx;  fprintf('  %4.2f',xx(i));  end;
fprintf('\nY (n=%2d):',ny);
for i=1:ny;  fprintf('  %4.2f',yy(i));  end;

disp(' ');   % For vertical spacing
SampleSizes_X_Y_Both = [ nx ny nn ]
%% fprintf('\n\nSample sizes:  %d  %d    Together: %d\n',nx,ny,nn);

%% To find joint ranks, we need a nnxd matrix with data in one column
%%   and group numbers in another column
xxcol = [ones(nx,1)  xx'];    %' X values with group number=1
yycol = [2*ones(ny,1)  yy'];  %' Y values with group number=2

% Data structure
dd = [xxcol;  yycol];         % All data in an nnx2 matrix
dd = [(1:nn)'  dd];           %' Prepend a row number

disp('Columns are Ordinal, Group number, and Values:');
%% Display dd (disp(dd); would suppress the  dd =   line)
dd

% Find midranks and tie-group sizes with tsize >= 2
zdat = dd(:,3);      % Column 3 is data values from both groups
% Call getmidranks() to find midranks and tie-group sizes
[midr tiegr] = getmidranks(zdat);
% Add midranks as a 4th column to dd
dd = [dd  midr];

disp('Tie-group sizes for tie-groups with more than one value:');
g=0;    % Find the number of tie-group sizes
% First end for if statement; second end for loop
for i=1:nn;  if tiegr(i)>1 g=i;  else break;  end;  end;
if g==0   Number_of_Tie_Groups = g    % Display
else disp(tiegr(1:g)');  end;    %' Display as row vector

disp('Columns are Ordinal, Group numbers, Values, and Midranks:');
disp('(The function  getmidranks()  was used to generate midranks.)');
disp('Note: Could use fprintf() to suppress unnecessary 0s');
dd   % Display

disp(' ');
disp('Data sorted by Z-value to show midranks and any tie groups');
disp('(Sorted on Z-value then Group number then Ordinal as tie-breakers)');
disp('Columns are Ordinal, Group number, Value, and Midrank');
disp('Using fprintf rather than matrix display for variety');
disp(' ');
ddx = sortrows(dd, [3 2 1]);
%% NOTE: To display by row with fprintf, must use ddx' (transpose=') not ddx
fprintf('   %2d.   %d   %4.2f   %5.1f\n', ddx');  %'

% All following calculations could be done using either dd or ddx  

% Calculate the Wilcoxon scores for both samples
% Column 2 has group numbers; column 4 has midranks
grp=dd(:,2);  midr=dd(:,4);
ObsWscoreX=0;  ObsWscoreY=0;
for i=1:nn;  if grp(i)==1 ObsWscoreX=ObsWscoreX+midr(i);
             else if grp(i)==2 ObsWscoreY=ObsWscoreY+midr(i);  end; end; end
% Expected values under random permutaions
WExpScX=nx*(nx+ny+1)/2;
WExpScY=ny*(nx+ny+1)/2;
% Base observed score on second (treatment) sample
Wdiff=ObsWscoreY-WExpScY;

disp('  ');
disp('WILCOXON RANK SUM SCORES:');
disp('(Since the treatment group (Y) is smaller, we show its score first)');
fprintf('Yscore:  Observed: %7.2f    Expected: %7.2f\n', ObsWscoreY,WExpScY);
fprintf('Xscore:  Observed: %7.2f    Expected: %7.2f\n', ObsWscoreX,WExpScX);
fprintf('Score-Expected in both cases:  %.2f\n', Wdiff);

disp('ASSUMING THAT THERE ARE NO TIES');
disp('(See WcxRankSum2.m for handling of ties)');

% Large-sample normal approximation
% Assuming NO TIE CORRECTION
Wvar = nx*ny*(nx+ny+1)/12;
Wstdev = sqrt(Wvar);
Zval = Wdiff/Wstdev;
Pval = normpval(Zval);

Wstd_Zval_Pval_NOTIECORR = [Wstdev  Zval  Pval]

% Carry out simulations to estimate true P-value

% Initialize random-number generator
wseed=sum(100*clock);

disp('SIMULATING THE TRUE Wilcoxon Rank-Sum P-VALUE:');
disp('  ');
disp('Using NSIMS random permutations:');

Nsims=10000

disp('INITIALIZING the random-number generator at');
Wseed=sum(100*clock())
rand('state',Wseed);

disp('ASSUMING ObsWscoreY is on the LOWER TAIL:');
disp('(That is, ObsWscoreY should be LESS than WExpScY)');
ObsWscoreY_ExpWscoreY = [ObsWscoreY  WExpScY]

% In the loop, kk is a random permutation of 1:nn
NLessEq=0;
for ns=1:Nsims;  kk=randperm(nn);
  % Note sum from Nx+1 to NN=Nx+Ny only (the Y sample)
  Yrand=0;  for i=nx+1:nn;  Yrand=Yrand+midr(kk(i));  end;
  if Yrand <= ObsWscoreY  NLessEq=NLessEq+1;  end;  end

phat=NLessEq/Nsims;     % one-sided P-value
SimPval = 2*phat;       % two-sided P-value

disp(' ');
disp('SIMULATED value for true 2-sided P-value:')

fprintf('NLessEq=%d  Nsims=%d   Pvalue=%7.5f (2-sided)\n',NLessEq,Nsims,SimPval);

disp(' ');
disp('95%% Confidence Interval for true P-value (with Estimate in middle):');
% phat, not SimPval, is distributed as Binom(N,p)/N 
pwid=1.960*sqrt(phat*(1-phat)/Nsims);
Plow=2*(phat-pwid);  Phigh=2*(phat+pwid);

Plow_PEst_Phigh = [Plow  SimPval  Phigh]

disp(' ');
disp('The Hodges-Lehmann estimator for theta for Y=X+theta')
disp('  is the median of the Mann-Whitney differences W_k=Y_i-X_j:')
% Collect the Mann-Whitney differences in a row vector 
ww=zeros(1,nx*ny+3);
r=0;
for i=1:nx;  for j=1:ny;  r=r+1;
   ww(r)=yy(j)-xx(i);  end;  end;
nmwd=r;   % Number of Mann-Whitney differences

% Display Walsh averages unsorted then sorted
nrows=ceil(nmwd/10);
fprintf ('\nMann-Whitney differences (n=%d, nmwd=%d, nrows=%d):\n', nn,nmwd,nrows);
for i=1:nrows;  r=(i-1)*10;  r2=r+10;
  if r2>nmwd r2=nmwd;  end;
  if r+1 > nmwd break;  end;  % end of if statement
  for j=r+1:r2;  fprintf('  %5.2f', ww(j));  end;
  fprintf('\n');  end;

fprintf ('\nSORTED M-W differences (n=%d, nmwd=%d, nrows=%d):\n', nn,nmwd,nrows);
ww=sort(ww(1:nmwd));
for i=1:nrows;  r=(i-1)*10;  r2=r+10;
  if r2>nmwd r2=nmwd;  end;
  if r < nmwd fprintf(' (%3d)',r+1);
  else break;  end;  % of if statement
  for j=r+1:r2;  fprintf('  %5.2f', ww(j));  end;
  fprintf('\n');  end;


% Display Zmean and MWMed
xmean=sum(xx)/nx;  ymean=sum(yy)/ny;
dmean=ymean-xmean;        % Sample Ybar-Xbar
MWMed = median(ww);
YXmeandiff_MWMedian = [dmean  MWMed]

fprintf ('\nNormal confidence intervals for theta (Mean diff=%g):\n',dmean);
xc=xx-xmean;  yc=yy-ymean;
svar = ((1/nx) + (1/ny))*(sum(xc.*xc) + sum(yc.*yc))/(nn-2);
sdev = sqrt(svar);
cwid = 1.960*sdev;
fprintf ('  (%6.3f, %6.3f)  95%% CI\n', dmean-cwid, dmean+cwid);
cwid = 2.576*sdev;
fprintf ('  (%6.3f, %6.3f)  95%% CI\n', dmean-cwid, dmean+cwid);

disp(' ');
disp('EXACT NONPARAMETRIC CIs for theta=med(Y)-med(X)');
disp('  based on the Wilcoxon rank-sum statistic:');

% From Table A.6 for n=5 m=10:
TabA6 = [ 61 0.004;  60 0.006;  59 0.010;
          58 0.014;  57 0.020;  56 0.028];
nTab=size(TabA6,1);  % Number of entries

maxScore=ny*(2*nx+ny+1)/2;
maxMWD=ny*nx;
nn2=ny*(ny+1)/2;
fprintf('\nTable A6 for n=%d m=%d  (MaxScore=%d)  has:\n',ny,nx,maxScore);
% Display TabA6
fprintf ('   P(W_Y ge %2d)=%5.3f\n', TabA6');  %'
fprintf ('\nWith corresponding Mann-Whitney scores:\n');
MWTab = [ TabA6(:,1)  TabA6(:,1)-nn2  TabA6(:,2)];
fprintf ('   P(W_Y ge %2d) = P(MW_Y ge %2d) = %5.3f\n', MWTab');  %'

disp(' ');
disp('Offsets in Mann-Whitney array with coverage probabilities:');
xhi=TabA6(:,1)-nn2;  xlo=nx*ny + 1 - xhi;
clevel=100*(1-2*TabA6(:,2));

ValTab = [ xlo  xhi  clevel ];
fprintf('  Lo=%4.1f   Hi=%4.1f   Cprob=%5.2f\n', ValTab');  %'

disp('  ');
disp('Exact symmetric confidence intervals for theta:');
CITab = zeros(nTab,2);
for i=1:nTab;
  CITab(i,1)=ww(xlo(i));  CITab(i,2)=ww(xhi(i));  end;
CITab = [CITab  clevel];

fprintf('  (%6.3f, %6.3f)   %4.1f%% CI\n', CITab');  %'

