function [pvalue,W,U,Wstar]=wilcoxranksum(X,Y,Ex)
% Wilcoxon rank-sum test
% [pvalue,W,U,Wstar]=wilcoxranksum(X,Y,Ex)
% Tests the null hypothesis that X & Y have the same pdf.
% Input: X,Y two samples,Ex~=0 indicates do an exact test.
% Output: pvalue: pvalue, 2-sided p value for large sample approximation N(0,1) distribution
%         W=Wilcoxon rank sum statistic
%         U=Mann-Whitney U statistic
%         Wstar=z value for asymptotic large sample approximation
%                Don't request Wstar for exact test.
% Calls Wilcoxrsexact
% Written by Eugene.Gallagher@umb.edu
% Revised 11/14/10
if nargout==4
    error('Wstar can''t be calculated for exact test')
    return
end
X=X(:);Y=Y(:);
n=length(X);
m=length(Y);
% Rank the X&Y values from smallest to largest, assigning average ranks to ties.
[T,R,ind]=ties([X;Y]);T=T';  % calls Gallagher's ties.m
% Find sum of ranks of the smaller sample;
if n<m;
    W=sum(R(1:n));
else
    W=sum(R(n+1:n+m));
    n=m;  % Expected value & variance equastions assume n is the size of the smaller group.
    m=length(X);
end
U=W-n*(n+1)/2; % Mann-Whitney U statistic
largesample=logical(1);
if nargin>2
    if Ex~=0
        largesample=logical(0);
    end
end
if nargin>2 & ~largesample
  ncomb=nchoosek(n+m,n);  
  if ncomb>1e6
      t=sprintf(...
          '%d combinations, T=%d min (1e6 combs take 1 min on p4)\n',...
          ncomb,round(ncomb/1e6)); 
      toomany=menu(t,'Stop','Continue');
      if toomany==1
         largesample=logical(1);fprintf('Large sample approximation for 2-tailed p\n');
      end
  end
  if ~largesample
      pexuptail=wilcoxrsexact(n,m,W,R);
      if pexuptail<=0.5
          pvalue=2*pexuptail;
      else
          pvalue=2*(1-pexuptail);
      end
  end
end
if largesample
    % Large sample approximation;% Hollander & Wolfe p. 108
    EoW=(n*(m+n+1))/2;
    % Calculate the variance of W, without ties and with ties.
    if isempty(T)  % Size of tied groups from ties.m
        VaroW=(m*n*(m+n+1))/12;
    else
        VaroW=(m*n)/12*(m+n+1-(sum((T-1).*T.*(T+1)))/((m+n)*(m+n-1)));
    end
    Wstar=(W-(n*(m+n+1)/2))/sqrt(VaroW);  % Without ties, tends to an asymptotic N(0,1) distribution.
    % Find the 2-tailedprobability of Wstar from the standard normal distributioin
    pvalue=erfc(abs(Wstar)/sqrt(2));
    % Note that the exact p values are tabulated, and an exact test, even in the presence of ties
    % can be performed, see pp. 113-116 in Hollander & Wolfe.
end

function [T,R,ind]=ties(A)
% format: [T,R,ind]=ties(A)
% a function to return a row vector of tied groups, T,
% Ranks R (including average ranks) and indices of tied elements
% needed to calculate variance of S using Kendall's 
% variance formula & Spearman's r.
% input: A is a row or column vector
% T: a row vector containing number of members of tied groups
% T=0 if there are no tied groups
% sum(T) is equal to the number of tied elements.
% each element of T equals the number in each tied group
% tied groups are sorted in ascending order.
% Examples: A=[1 2 3];[T,R,ind]=ties(A)=>   T=0,R=[1 2 3],ind=[]
%           A=[1 2 3 1];           T=2,R=[1.5 3 4 1.5],ind=[1 4]
%           A=[2 1 2 3 1 2];       T=[2 3],R=[4 1.5 4 6 1.5 4],
%                                  ind=[5 2 3 1 6]
%           A=[2 1 2 3 3 1 2];     T=[2 3 2],R=[4 1.5 4 6.5 6.5 1.5 4]
%                                  ind=[6 2 3 1 7 4 5]
% R (Row vec)=numerical rankings of A with ave. ranks for ties
% ind: indices of tied elements, sorted by rank; sorted tied elements=A(ind);
% ties.m is used in Kendall.m as T=ties(A), and Spear.m
% written by E. Gallagher, Environmental Sciences Program 
% UMASS/Boston, Email: Eugene.Gallagher@umb.edu
% written: 6/16/93, revised 6/17/93
[r,c]=size(A);
if r>c
     A=A';                   % change to row vector
end
[Asort,k]=sort(A);
iota=1:length(A);iota=iota';
R(k)=iota;
index=[k' iota];
ind=[];
CDA=[~diff(Asort) 0];
min1=min(find(CDA==1));
if isempty(min1)
     T=0;
     return
end
i=0;
[rw,cl]=size(CDA);
T=zeros(size(rw,cl));
while ~isempty(min1)
     min0=min(find(CDA==0));
     if min0<min1
           CDA(min0:min1-1)=[];
           index(min0:min1-1,:)=[];
     else
       i=i+1;
       T(i)=min0-min1+1;
       CDA(min1:min0)=[];
       ind=[ind index(min1:min0,1)'];
       R(1,index(min1:min0))=ones(1,T(i))*sum(index(min1:min0,2))/T(i);
       index(min1:min0,:)=[];
     end
     min1=min(find(CDA==1));
end
T(find(T==0))=[];