function [tau,L] = kendall_tau(X)

%KENDALL_TAU computes Kendall's tau statistic for two distributions.

%INPUTS:  X = nx2 data matrix

%OUTPUT: tau = Kendall's tau

n = size(X,1);

%*******************************

%First compute ranks

%*******************************

x1 = X(:,1);

x2 = X(:,2);

L1 = [[1:n]',x1];

L1 = sortrows(L1,2);

R1 = zeros(n,1);

L2 = [[1:n]',x2];

L2 = sortrows(L2,2);

R2 = zeros(n,1);

%Now convert x1 to proper ranks

rows = L1(:,1);

R = [1:n]';

L = L1(:,2);

i = 1;

D1 = n*(n-1)/2;

while i < n
   
   ctr = i;
   
   j = 1;
   
   while (i + j <= n)&(L(i+j) == L(i))
      
      ctr = ctr + (i + j);
      
      j = j + 1;
      
   end
   
   if j > 1
      
      R(i:i+j-1) = ctr/j;  %Revise ranks
      
      D1 = D1 - j*(j-1)/2; %Revise D1
      
      i = i + j;
      
   else
      
      i = i + 1;
      
   end
   
end

R1(rows) = R;


%Next convert x2 to proper ranks

rows = L2(:,1);

R = [1:n]';

L = L2(:,2);

i = 1;

D2 = n*(n-1)/2;

while i < n
   
   ctr = i;
   
   j = 1;
   
   while (i + j <= n)&(L(i+j) == L(i))
      
      ctr = ctr + (i + j);
      
      j = j + 1;
      
   end
   
   if j > 1
      
      R(i:i+j-1) = ctr/j;  %Revise ranks
      
      D2 = D2 - j*(j-1)/2; %Revise D2
      
      i = i + j;
      
   else
      
      i = i + 1;
      
   end
   
end

R2(rows) = R;


%*******************************

%Now compute Kendall's Tau

%*******************************

S = 0;

for i = 1:n-1
   
   x1 = R1(i);
   
   y1 = R2(i);
   
   for j = i+1:n
      
      x2 = R1(j);
      
      y2 = R2(j);
      
      S = S + sign(x1 - x2)*sign(y1 - y2);
      
   end
   
end