function [p, h, stats] = signtest(x,y,alpha)
%SIGNTEST Sign test.
%   P = SIGNTEST(X,Y,ALPHA) returns the significance a test of the null
%   hypothesis that the median difference between two matched samples,
%   X and Y, is zero.  X and Y must be vectors of equal length. Y may
%   also be a scalar.  In this case, SIGNTEST computes the probability
%   that the median of X is different from the constant, Y.
%
%   ALPHA is the desired level of significance.  It must be a
%   scalar between zero and one, and its default value is 0.05.
%
%   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] = SIGNTEST(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] = SIGNTEST(X,Y,ALPHA) also returns STATS, 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
%   Copyright 1993-2000 The MathWorks, Inc. 
% $Revision: 1.10 $

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

[rowx, colx] = size(x);
[rowy, coly] = size(y);
if prod(size(y)) == 1,
   rowy = max(size(x));
   y = y(ones(rowy,1));
end

if min(rowx, colx) ~= 1 | min(rowy,coly) ~= 1,
   error('SIGNTEST 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('SIGNTEST requires the data vectors to have the same number of elements.');
end

n = rowx;
diffxy = x - y;
nodiff = find(diffxy == 0);
diffxy(nodiff) = [];
n = n - length(nodiff);

if n == 0 % this means the two vectors are identical
    p = 1;
    if nargout == 2,
	h = (p<alpha);
    end   
    return
end

npos = length(find(diffxy>0));
nneg = n-npos;
sgn = min(nneg,npos);

if n < 100
   p = min(1, 2*binocdf(sgn,n,0.5));  % p>1 means center value double-counted
else
   % Do a continuity correction, keeping in mind the right direction
   z = (npos-nneg - sign(npos-nneg))/sqrt(n);
   p = 2*normcdf(-abs(z),0,1);
   if (nargout > 2)
      stats.zval = z;
   end
end

if nargout > 1
   h = (p<=alpha);
   if (nargout > 2)
      stats.sign = sgn;
   end
end   
