function [p, h, stats] = signrank(x,y,alpha)
%SIGNRANK Wilcoxon signed rank test of equality of medians.
%   P = SIGNRANK(X,Y,ALPHA) returns the significance for a test of the
%   null hypothesis that the median difference between two matched
%   samples, X and Y, is zero.
%
%   ALPHA is the desired level of significance, and must be a
%   scalar between 0 and 1.  Its default value is 0.05.
%
%   If Y is a scalar, it is extended to the length of X.  This
%   facilitates comparing the median of one sample to a constant.
%
%   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.
%
%   [P,H] = SIGNRANK(X,Y,ALPHA) also returns H, the result of the
%   hypothesis test.  H is 0 if the medians of X and Y are not
%   significantly different, and 1 if they are significantly
%   different.
%
%   [P,H,STATS] = SIGNRANK(X,Y,ALPHA) also returns a STATS is a
%   structure with one or two fields.  The field 'signedrank' contains
%   the value of the signed rank statistic.  If the sample size is
%   large, then P is calculated using a normal approximation and the
%   field 'zval' contains the value of the normal (Z) statistic.

%   B.A. Jones 11-18-96, ZP You 9-4-98
%   Copyright 1993-2000 The MathWorks, Inc. 
% $Revision: 1.11 $

if nargin < 3
   alpha = 0.05;
end
if (length(alpha)>1)
   error('SIGNRANK requires a scalar ALPHA value.');
end
if (isnan(alpha) | (alpha <= 0) | (alpha >= 1))
   error('SIGNRANK requires 0 < ALPHA < 1.');
end

[rowx, colx] = size(x);
if (length(y) == 1)
   y = repmat(y, rowx, colx);
end
[rowy, coly] = size(y);

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,
   error('SIGNRANK requires the data vectors to have the same number of elements.');
end

diffxy = x - y;

diffxy(isnan(diffxy)) = [];
if (length(diffxy)==0), error('No data remaining after removal of NaNs.'), end

nodiff = find(diffxy == 0);
diffxy(nodiff) = [];
n = length(diffxy);

if (n == 0)         % degenerate case, all ties
   p = 1;
   if (nargout > 1)
      h = 0;
      if (nargout > 2)
         stats.signedrank = 0;
      end
   end
   return
end

% Find negative differences and ranks of absolute differences
neg = find(diffxy<0);
[tierank, tieadj] = tiedrank(abs(diffxy));

% Compute signed rank statistic (most extreme version)
w = sum(tierank(neg));
w = min(w, n*(n+1)/2-w);

% Compute significance exactly or with a normal approximation
if n > 15,
   z = (w-n*(n+1)/4) / sqrt((n*(n+1)*(2*n+1) - tieadj)/24);
   p = 2*normcdf(z,0,1);
   if (nargout > 2)
      stats.zval = z;
   end
else
   allposs = (ff2n(n))';
   idx = (1:n)';
   idx = idx(:,ones(2.^n,1));
   pranks = sum(allposs.*idx);
   tail = 2*length(find(pranks <= w)); % two side.

   % Avoid p>1 if w is in the middle and is double-counted
   p = min(1, tail./(2.^n));
end

if nargout > 1,
   h = (p<=alpha);
   if (nargout > 2)
      stats.signedrank = w;
   end
end
