%
% Large-sample approximation and simulation of true P-values
%   for Wilcoxon rank-sum test
% Use  getmidranks  function to find midranks
% Use normal approximation for order-statistic offsets
%   for nonparametric confidence interval

echo off;

disp('Wilcoxon Rank-Sum test for two samples:');

xx = [23  21  22  14  23  19  15  15  16];
yy = [27  27  22  21  27  28  19  20  31  24  21  23  23  21];

nx=size(xx,2);    % Sample sizes
ny=size(yy,2);

% Display the two samples
fprintf('\nX:');
for i=1:nx;  fprintf('  %d',xx(i));  end;
fprintf('\nY:');
for i=1:ny;  fprintf('  %d',yy(i));  end;

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

%% Adjoin group numbers to treatment-group data columns
%%   to store data from both samples in a single
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 (could also say  disp(dd);)
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('The function  getmidranks()  was used to generate midranks');
disp('  and tie-group sizes');
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;

Number_TieGroups_of_Size_ge_2 = g     % Display
if g~=0 
   disp('TieGroup Sizes of Size ge 2:');  disp('  ');
   disp(tiegr(1:g)');  end;    %' Display as row vector

disp('Columns are Ordinal, Group numbers, Values, and Midranks:');
disp('Note: Directly displaying matrix dd would have all columns');
disp('  to 4 places beyond decimal point because of midranks.');
disp('Instead, we use fprintf() on the rows of dd so that we can');
disp('  determine the format.');
disp(' ');
% Note dd' (dd transpose) in syntax for fprintf(), NOT dd    %'
% Since midranks can be fractions, the usual MATLAB display would give
%   all numbers to 4 decimal places
fprintf('   %2d.   %d   %2d   %5.1f\n', dd');  %'  

disp(' ');
disp('Data sorted by Z-value to show 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(' ');
ddx = sortrows(dd, [3 2 1]);
fprintf('   %2d.   %d   %2d   %5.1f\n', 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 first sample
Wdiff=ObsWscoreX-WExpScX;

disp('  ');
disp('WILCOXON RANK SUM SCORES:');
fprintf('Xscore:  Observed: %7.2f    Expected: %7.2f\n', ObsWscoreX,WExpScX);
fprintf('Yscore:  Observed: %7.2f    Expected: %7.2f\n', ObsWscoreY,WExpScY);
fprintf('Difference in both cases:  %.2f\n', Wdiff);

% Large-sample normal approximation with NO TIE CORRECTION
WvarNoTC = nx*ny*(nx+ny+1)/12;
WstdevNoTC = sqrt(WvarNoTC);
ZvalNoTC = Wdiff/WstdevNoTC;
PvalNoTC = normpval(ZvalNoTC);

Wstd_Zval_Pval_NOTIECORR = [WstdevNoTC  ZvalNoTC  PvalNoTC]

%% Compute tie-sum correction from tie-group sizes
Tiesum=0;
for i=1:nn;  tt=tiegr(i);
  if tt>1 Tiesum=Tiesum+(tt-1)*tt*(tt+1);
  else break;  end;  end;
TieCorr = Tiesum/((nn-1)*nn*(nn+1));

Tiesum_TieCorr = [Tiesum  TieCorr]

% Tie correction to normal P-value
WvarTC=WvarNoTC*(1-TieCorr);
WstdevTC = sqrt(WvarTC);
ZvalTC = Wdiff/WstdevTC;
PvalTC = normpval(ZvalTC);

Wstd_Zval_Pval_TIECORR = [WstdevTC  ZvalTC  PvalTC]

% Carry out simulations to estimate true P-value

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

disp('SIMULATING THE TRUE Wilcoxon Signed-Rank 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 ObsWscoreX is on the LOWER TAIL:');
disp('(That is, ObsWscoreX should be LESS than WExpScX)');
ObsWscoreX_ExpWscoreX = [ObsWscoreX  WExpScX]

% In the loop, kk is a random permutation of 1:nn
NLessEq=0;
for ns=1:Nsims;  kk=randperm(nn);
  % Note sum from 1 to Nx only
  Wrand=0;  for i=1:nx;  Wrand=Wrand+midr(kk(i));  end;
  if Wrand <= ObsWscoreX  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('  %4.0f', 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('  %4.0f', 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]


disp('APPROXIMATE EXACT CONFIDENCE INTERVALS');
disp('  based on the Wilcoxon rank-sum statistic:');

fprintf('\nEstimator: %g\n', MWMed);

% In contrast, normal confidence intervals:
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 ('  (%5.2f, %5.2f)  95%% CI\n', dmean-cwid, dmean+cwid);
cwid = 2.576*sdev;
fprintf ('  (%5.2f, %5.2f)  99%% CI\n', dmean-cwid, dmean+cwid);

disp(' ');
disp('NONPARAMETRIC approximate exact CIs for theta');
disp('  based on the Wilcoxon rank-sum statistic:');
% Offsets in ww() array for nx=9, ny=14:
cmean=nx*ny/2;
cstd=sqrt(nx*ny*(nn+1)/12);
% Approximate 95%:
r1 = floor(cmean - 1.960*cstd);
r2 = nx*ny+1 - r1;
fprintf ('  (%2.0f, %2.0f)  Offsets r1=%2d  r2=%3d  95%% CI\n', ww(r1),ww(r2),r1,r2);
% Approximate 95%:
r1 = floor(cmean - 2.576*cstd);
r2 = nx*ny+1 - r1;
fprintf ('  (%2.0f, %2.0f)  Offsets r1=%2d  r2=%3d  99%% CI\n', ww(r1),ww(r2),r1,r2);


