function [p, h] = signrankp(x,y,alpha)
%SIGNRANK Wilcoxon signed rank test of equality of means for
%comparing samples of unequal size.
%
%   P = SIGNRANK(X,Y,ALPHA) returns the significance probability
%   that the means of two samples, X and Y are equal.
%   X and Y need not be vectors of equal length. ALPHA is the desired
%   level of significance. and must be a scalar between
%   zero and one.
%
%   [P, H] = SIGNRANK(X,Y,ALPHA) also returns the result of the 
%   hypothesis test, H. H is zero if the difference in means of
%   X and Y is not significantly different from zero. H is one if
%   the two means are significantly different. 
%
%   P is the probability of observing a result equally or more 
%   extreme than the one using the data (X and Y) if the null  
%   hypothesis is true. If P is near zero, this casts doubt on
%   this hypothesis.
%
%   Currently works for sample sizes > 25.
%
% see also SIGNRANK, RANKSUM
%
% Copyright(c): PNath@London.edu 12-Mar-2001
%

if nargin < 3
   alpha = 0.05;
end

[rowx, colx] = size(x);
[rowy, coly] = size(y);

if min(rowx,rowy) < 25
   error('SignRankP currently only works for sample sizes > 25');
end


if min(rowx, colx) ~= 1 | min(rowy,coly) ~= 1,
   error('SIGNRANK requires vector rather than matrix data.');
end 
if rowx == 1
   rowx = colx;
   x = x';
end
if rowy == 1,
   rowy = coly;
   y = y';
end
   
if rowx == rowy,
   [p, h] = signrank(x,y,alpha);
   return
end

CombinedSample = [x;y];
CombinedSample(:,2) = Order(CombinedSample,2);   % returns rank adjusted for non-uniqueness in 2nd col
RankedX = CombinedSample(1:rowx, 1:2);
RankedY = CombinedSample(rowx+1:rowx+rowy, 1:2);
if rowx < rowy
   T1 = sum(RankedX(:,2));
   T2 = rowx*(rowx + rowy + 1) - T1;
   T = min(T1,T2);
   MuT = rowx*(rowx+1)/4;
   SigmaT = (rowx*(rowx+1)*(2*rowx+1)/24)^.5;
   ZStat = (T-MuT)/SigmaT;
else
   T1 = sum(RankedY(:,2));
   T2 = rowy*(rowx + rowy + 1) - T1;
	T = min(T1,T2);
   MuT = rowy*(rowy+1)/4;
   SigmaT = (rowy*(rowy+1)*(2*rowy+1)/24)^.5;
   ZStat = (T-MuT)/SigmaT;
end

p = 1-normcdf(abs(ZStat));

 


if 1==2     %this part from singrank is not used
diffxy = x - y;
nodiff = find(diffxy == 0);
diffxy(nodiff) = [];
n = length(diffxy);
[sd, rowidx] = sort(abs(diffxy));
neg = find(diffxy<0);

invr(rowidx) = 1:n; % invr is the inverse of rowidx.
w = sum(invr(neg));
w = min(w, n*(n+1)/2-w);

if n > 15,
   z = (w-n*(n+1)/4)/sqrt(n*(n+1)*(2*n+1)/24);
   p = 2*normcdf(z,0,1);
else
   allposs = (ff2n(n))';
   idx = (1:n)';
   idx = idx(:,ones(2.^n,1));
   pranks = sum(allposs.*idx);
   tail = 2*length(find(pranks < w))+length(find(pranks == w)); % two side.
   p = tail./(2.^n);
end
end


if nargout == 2,
   h = (p<alpha);
end