function  [multcomp]=friedmanpairwise( data, alpha  )
% ------------------------------------------------------------------------------------
%   FRIEDMAN's  PAIRWISE COMPARISONS 
%   Input:  data  - data matrix in which rows are blocks and columns are treatments. 
%           alpha - level, default 5%
%
%   Output:
% CI's for sum rank differences. If the CI contains 0
% there is no significant difference between the 
% corresponding treatments. This code follows text 
% and Conover's Practical NP Statistics (Wiley, 1999)
% page 371.
%                      
%   Comes with ranks.m (ranking procedure) included
%   Ver 0.3 July 2014
%-------------------------------------------------------------------------------------
if nargin < 2
    alpha = 0.05;
end
[b, k] = size(data);
ranksj = sum(ranks(data',0)');     % Here we find the R_j's
a = sum(sum( ranks(data',0).^2 ));
df = (b-1)*(k-1);
pt = tinv( 1-alpha/2 , df );
multcomp=[];
displace= pt * sqrt( 2*(b*a - sum(ranksj.^2))/df);
for i=1:k-1
    for j=i+1:k
     difij = ranksj(i) - ranksj(j);
multcomp=[multcomp; 
    i j difij - displace  difij  difij + displace];
    end
end
%

% [b, k] = size(data);
% ranksj = sum(ranks(data',0)');
% a = sum(sum( ranks(data',0).^2 ));
% c = b*k*(k+1)^2/4;
% %--------------------------------------------------------------------
% S = (k-1)* ( sum( ranksj.^2 ) - b*c )/(a -c);
% F = (b-1)*S/(b*(k-1)-S);



%----------------------------------------------------------------------
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