function  [W, Z, p] = wsurt( data1, data2, alt )
% ------------------------------------------------------------------------------------
%   WILCOXON SUM RANK TEST
%   Input:  data1, data2 - first and second sample
%           alt - code for alternative hypothesis;
%               -1  mu1<m2; 0 mu1 ~= m2; and 1 mu1>mu2
%   Output: W - sum of the ranks for the first sample. If
%               there is no ties, the standardization by EW = ... and
%               Var W =  ... , allows using standard normal quantiles as long as
%               sample sizes are larger than 15-20.
%           Z - standardized R but adjusted for the ties
%           p - p-value for testing equality of distributions (equality of locations)
%               against the alternative specified by input alt
%   Example of use:
%   >  dat1=[1 3 2 4 3 5 5 4 2 3 4 3 1 7 6 6 5 4 5 8 7 3 3 4];
%   >  dat2=[2 5 4 3 4 3 2 2 1 2 3 2 3 4 3 2 3 4 4 3 5];
%   >  [sumranks1, tstat, pval] = wsurt(dat1, dat2, 1)
%
%   Comes with function ranks.m (ranking procedure)
%-------------------------------------------------------------------------------------
data1 = data1(:)' ;                   % convert sample 1 to a row vector
  n1 = length( data1 );               % n1  is the size of the first sample, data1
data2 = data2(:)' ;                   % convert sample 2 to a row vector
  n2 = length( data2 );               % n2 is the size of the second sample, data2
  n=n1+ n2;                           % n is the total sample size
mergeboth = [ data1 data2 ];
rankall = ranks( mergeboth );       % ranks of merged observations
  W2 = sum( ( rankall.^2 ) );       % sum of all ranks squared needed to
% make adjustment for the ties; if no ties are present, this sum is equal
% to the sum of squares of the first n integers: n(n+1)(2 n+1)/6.
ranks_data1 = rankall( :, 1:n1);    % ranks of the first sample only
  W = sum( ranks_data1 );             % statistics for WSuRT, sum of the ranks for sample 1
%--------------------------------------------------------------------------------------
Z = (W - n1 * (n+1)/2 )/sqrt( n1*n2*W2/(n * (n-1)) - n1*n2*(n+1)^2/(4*(n-1)));
% Z is approximately standard normal and approximation is quite good if n1,n2 > 15.
% Since W ranges over integers and half integers, a continuity correction, cc, may be
%  used for improving accuracy of p-values.
cc = 0.25/sqrt( n1*n2*W2/(n * (n-1)) - n1*n2*(n+1)^2/(4*(n-1)));
%-------------------------- alternatives ---------------------------------------
% alt == 0 for two sided;  alt == -1 for mu1 < mu2; alt == 1 for mu1 > mu2.
if      alt == 0
            p = 2*normcdf(- abs( Z ) - cc);
elseif  alt == -1
            p = normcdf( Z - cc );
elseif  alt == 1
            p = normcdf( - Z - cc );
else
        error('Wrong alternative code. Input "alt" should be either 0, -1, or 1.')
end
%-------------------------------------------------------------------------------------
% Copyright (C) 2002 by GaTech
% Ver 0.1   9/28/2002
%----------------------------------------------------------------------------------
function  r = ranks(data, glob)
%--------------------------------------------------------------------
% RANK TRANSFORMATION of data on input
%    Input: data - data to be ranked, data is matrix
%           glob - how to rank data matrix. glob=1 ranks amm matrix components (global 
%                  ranking; glob ne 1 - rank by columns of data. Default glob=1.
%    Output:   r - rank-transformed data
% Finds ranks of a vector of observations and adjusts for ties (by averaging)
% 
% Comment. If ties are no present (ties in ranking are matching observations), 
% ranks can be found simply as the indices in sorting of the indices of the
% sorted data (sentence is OK), i.e.,  [a,b]=sort(data); [c,d]=sort(b);
% and the ranks are in 'd'. If there are ties, they are ranked from left to
% right and also from right to left, and then averaged!
%---------------------------------------------------------------------
if nargin < 2
    glob = 1;
end
    shape = size(data);
if glob == 1    
     data=data(:);
 end     
% Ties ranked from UptoDown
  [ irrelevant , indud ]  =  sort(data);
  [ irrelevant , rUD ]    =  sort(indud);
% Ties ranked from RtoL
  [ irrelevant , inddu ]  =  sort(flipud(data));
  [ irrelevant , rDU ]    =  sort(inddu);
% Averages ranks of ties, keeping ranks of no-tie-observations the same
r = (rUD + flipud(rDU))./2;
r = reshape(r,shape); 
%---------------------------------------------------------------------
% Copyright (C) 2002 by GaTech
% Ver 0.1   9/28/2002