function [pvalue,W]=wilcoxsignrank(X,Y)
% Wilcoxon signed rank test
% [pvalue,W]=wilcoxsignrank(X,Y)
% Tests the null hypothesis that X & Y have the same pdf.
% Input: X,Y two paired samples
% Output: pvalue:2-sided p value for large sample approximation N(0,1) distribution
%         W=Wilcoxon signed rank statistic, sum of the positive ranks.
% Calls Gallagher's ties.m (included below) to find the pattern of ties.
% Reference: Hollander & Wolfe (1999) Nonparametric statistical methods
%            Larsen & Marx (2001), page 692-695
% Written by E. Gallagher;
X=X(:);Y=Y(:);
n=length(X);m=length(Y);
if n~=m
    error('X & Y must be the same length')
end
% Rank the abs(X-Y) values from smallest to largest, assigning average ranks to ties.
Diff=X-Y;Psi=Diff>0;
[T,R,ind]=ties(abs(Diff));T=T';R=R';  % calls Gallagher's ties.m
W=sum(Psi.*R);   % Hollander & Wolfe, Equation 3.3.

% Large sample approximation;% Hollander & Wolfe p. 108
EoW=(n*(n+1))/4;
% Calculate the variance of W, without ties and with ties.
if isempty(T)  % Size of tied groups from ties.m
    VaroW=(n*(n+1)*(2*n+1))/24;
else
    VaroW=( n*(n+1)*(2*n+1) -0.5*(sum((T-1).*T.*(T+1)))    ) / 24;
end
Tstar=(W-EoW)/sqrt(VaroW);  % Without ties, tends to an asymptotic N(0,1) distribution.
% Find the 2-tailed probability of Wstar from the standard normal distribution
pvalue=erfc(abs(Tstar)/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.
