function  [table]= kwpairwise(data, belong, alpha)
% ------------------------------------------------------------------------------------
%   KRUSKAL-WALLIS PAIRWISE COMPARISONS
% Ver 0.1   4/5/2014
% Gives CI's for differences in average ranks between treatments i and j.
% If the CI contains 0, there is no significant difference between the
% treatments, if the interval for the difference is all positive
%(negative) then the  treatment i is superior (inferior) to treatment j.
%  >> data  =  [ 1 3 4  3 4 5  4 4 4 6 5 ];
%  >> belong = [ 1 1 1  2 2 2  3 3 3 3 3 ]; 
%  >> kwpairwise(data, belong)
%    1.0000    2.0000   -8.0097   -2.8333    2.3431
%    1.0000    3.0000   -9.1632   -4.5333    0.0966
%    2.0000    3.0000   -6.3299   -1.7000    2.9299
% *** 
%   Comes with functions ranks.m, and different.m  


if nargin<3
    alpha=0.05;
end
data = data(:)';      
  n = length( data );   
belong = belong(:)' ;      
  N = length( belong ); 
  if n ~= N 
      error('DATA and BELONG should be paired.')
  end
  Ra = ranks(data, 1);
  di = different(belong);
  k = length(di);
  ni=[];
 for i = 1:k
      ni= [ni sum( belong==di(i) ) ];
 end
  aver1  =[];   aver2 = [];
  for i = 1:k
      aver1= [aver1   sum(Ra(belong==di(i)))/ni(i)];
      aver2= [aver2   sum(Ra(belong==di(i)))^2/ni(i)];
  end
  S2= 1/(n-1) * (sum(Ra.^2) - 1/4*n*(n+1)^2);
Ti = sum(aver2);
H = 1/S2 * ( Ti -  1/4 * n*(n+1)^2);

table=[];
for (i = 1:k-1)
    for (j=i+1:k)
        ss2=S2*(n-1-H)/(n-k)*(1/ni(i)+ 1/ni(j));
        table=[table;  i j aver1(i) - aver1(j) - tinv(1-alpha/2,n-k)*sqrt(ss2) ...
         aver1(i) - aver1(j)    aver1(i) - aver1(j) + tinv(1-alpha/2,n-k)*sqrt(ss2) ];
    end
end
%-------------------------------------------------------------------------------------
% Copyright (C) 2002 by GaTech
% Ver 0.1   9/28/2002
% Ver 0.2   9/19/2003
%-------------------------------------------------------------------------------------
 function a = different(x)
 %
 % returns only different elements of x
 %
  if any (isnan (x))
    a = NaN;
  else
  a = [];
  x = x(:);
  s = sort(x);
  a = [s([1; find(s(2:end) > s(1:end-1)) + 1]); a];
  end
  %-------------
function  r = ranks(data, glob)
%--------------------------------------------------------------------
% 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