function  [T, p] = wsirtexa(data1, data2, alt )
% ------------------------------------------------------------------------------------
%   WILCOXON SIGNED RANK TEST EXACT
%   Input:  data1, data2 - first and second sample
%           data2 can be a scalar if the test is one sample
%           alt - code for alternative hypothesis;
%               -1  for mu1<mu2; 0 for mu1 ~= mu2; and 1 for mu1>mu2
%   Output: W - sum of all signed ranks
%           Z - standardized R but adjusted for the ties
%           p - p-value for testing the null hypothesis
%               against the alternative specified by the input alt
%   Examples of use:
%   >  dat1=[1 3 2 4 3 5 5 4 2 3 4 3 1 7 6 6 5 4 5 8 7];
%   >  dat2=[2 5 4 3 4 3 2 2 1 2 3 2 3 4 3 2 3 4 4 3 5];
%   >  [T, pval] = wsirtexa(dat1, dat2, 1)
%
 % When using wsirt as one sample test for mu=mu0 vs ~=
%  >  mu0 = 3.5;
%  >  dat1=[1 3 2 4 3 5 5 4 2 3 4 3 1 7 6 6 5 4 5 8 7];
%  > [T, pval] = wsirtexa(dat1,mu0, 0)
%  
%   Comes with function ranks.m (ranking procedure)

if nargin < 3
    alt=0
end
%-------------------------------------------------------------------------------------
data1 = data1(:)' ;       % convert sample 1 to a row vector
data2 = data2(:)' ;       % convert sample 2 to a row vector
if length(data2 == 1)  %for one sample test H0: mu=data2
      data2 = data2 .* ones(size(data1));
      %if data2 is constant repeat it to make vector data2
end
if length(data1) ~= length(data2);
    error('Sample sizes should coincide')
end
difs = data1 - data2;
difs = difs( difs ~= 0); %exclude 0 differences from consideration
n= length(difs);  nn = n*(n+1)/2;
rank_all = ranks(abs(difs));
   T = sum( rank_all .* (difs > 0));
%-------------------------- alternatives ---------------------------------------
% alt == 0 for two sided;  alt == -1 for mu1 < mu2; alt == 1 for mu1 > mu2.
if      alt == 0
            p  = 2* wsirtcdf(min(T, nn-T),n);
elseif  alt == -1
            p  = wsirtcdf(T,n);
elseif  alt == 1
            p = 1 - wsirtcdf(T,n);
else
        error('Wrong alternative code. Input "alt" should be either 0, -1, or 1.')
end
%-------------------------------------------------------------------------------------
% Copyright (C) 2002 by GaTech
% Ver 0.1   9/28/2002
function  r = ranks(data, glob)
%--------------------------------------------------------------------
% RANK TRANSFORMATION of data on input
%    Input: data - data to be ranked, data is matrix
%           glob - how to rank data matrix. glob=1 ranks amm matrix components (global 
%                  ranking; glob ne 1 - rank by columns of data. Default glob=1.
%    Output:   r - rank-transformed data
% Finds ranks of a vector of observations and adjusts for ties (by averaging)
% 
% Comment. If ties are no present (ties in ranking are matching observations), 
% ranks can be found simply as the indices in sorting of the indices of the
% sorted data (sentence is OK), i.e.,  [a,b]=sort(data); [c,d]=sort(b);
% and the ranks are in 'd'. If there are ties, they are ranked from left to
% right and also from right to left, and then averaged!
%---------------------------------------------------------------------
if nargin < 2
    glob = 1;
end
    shape = size(data);
if glob == 1    
     data=data(:);
 end     
% Ties ranked from UptoDown
  [ irrelevant , indud ]  =  sort(data);
  [ irrelevant , rUD ]    =  sort(indud);
% Ties ranked from RtoL
  [ irrelevant , inddu ]  =  sort(flipud(data));
  [ irrelevant , rDU ]    =  sort(inddu);
% Averages ranks of ties, keeping ranks of no-tie-observations the same
r = (rUD + flipud(rDU))./2;
r = reshape(r,shape); 
%---------------------------------------------------------------------
% Copyright (C) 2002 by GaTech
% Ver 0.1   9/28/2002


function [prob]=wsirtcdf(T, n)
% 
if T < 0  
    prob=0;
end

 if T >  n*(n+1)/2 
    prob = 1;
 end
 
 if T >= 0 && T <=  n*(n+1)/2 
 if n < 100 
    syms t
%probability generating function
     a=2^(-n)* prod( 1+t.^(1:n));
     co = coeffs(expand(a));
prob = double(sum(co(1:T+1))); %exact...+1 added since co starts with 0, not 1
else  
prob = normcdf(T+1+0.5, n*(n+1)/4,  sqrt(n*(n+1)*(2*n+1)/24)); %normal approximation
 end
 end
 
    
    
%---------------------------------------------------------------------
% Copyright (C) 2002 by GaTech
% Ver 0.1   9/28/2002
% Ver 1.0   10/3/2015


