% An experiment to determine the effect of digital massage and hyperosmolarity 
% on vitreous weight  was conducted by Galin et al. (1971). The results had 
% implications in  the context of cataract surgery. The authors measured
% the vitreous body weight for each eye of 15 New Zealand albino rabbits.
% One eye was treated (y), while the other served as a control (x).
% The measurements  are provided next:
% 
%           Galin, M. A., Robbins, R., and Obstbaum, S. (1971).
%           Prevention of vitreous loss.  Brit. J. Ophthal., 55, 533--537.

x=[1848 1532 1460 1947 1810 1718 1686 1617 ...
               1724 1873 1928 2226 1708 1605 1822];
y=[1738 1440 1388 1756 1692 1629 1583 1499 ...
               1596 1794 1785 2044 1602 1491 1702];

n=15

%nc-number of concordant pairs, 
%nd - discordant
%nt - ties
nc = 0; nd = 0; nt=0;
for i = 1:n-1
    for j = i+1:n
      if           sign(x(i) - x(j)) * sign(y(i) - y(j))  > 0     nc=nc+1;
      elseif    sign(x(i) - x(j)) * sign(y(i) - y(j))  < 0    nd=nd+1;  
      else   nt=nt+1;
      end
    end
end  
tauhat = (nc-nd)/(nc+nd)                       %0.9048
% OR 
%tau =corr(x', y', 'type','Kendall')  %0.9048
nc %100
nd %5
nt  %0
[u1 v1]=sort(x);
%[u2 v2]=sort(y);
v2=ranks(y);
s=0
z=v2(v1)
 for  i = 1:n-1 
     for j=i+1:n
         s=s+ (z(j) < z(i) );
     end
 end
  arara=1- 4*s/(n*(n-1))
  
%Sample variance of tau
tauvar = 2*(2*n+5)/(9*n*(n-1))     %0.0370
%tau is approx normal N(tau, 2*(2*n+5)/(9*n*(n-1))  )
%95% CI for population tau
[max(0,tauhat - norminv(0.975)*sqrt(tauvar)) ...
                  min(1,  tauhat + norminv(0.975)*sqrt(tauvar)) ]
%   [0.5276    1.0000]
S=nc-nd   %95
% p-value in testing H0: tau=0 vs H_1: tau != 0.
pval=2*(1-normcdf( abs(tauhat), 0, sqrt(tauvar))  )    %2.5853e-06
%or equivalently using Stau
pval1 = 2 *  (1-normcdf(3*abs(S)*sqrt(2)/sqrt(n*(n-1)*(2*n + 5))) )
 %2.5853e-06

 %Corrected two sided p-value (Stau is integer valued)
pval2 = 2 * min(normcdf(3*(S +1/2)*sqrt(2)/sqrt(n*(n-1)*(2*n + 5))) ,...   
             1-normcdf(3*(S -1/2)*sqrt(2)/sqrt(n*(n-1)*(2*n + 5))  ) )
% 2.9176e-06


%%
%Pearson corr coefficient
r=corr(x', y')  %0.9907
% A rough estimation of Pearson r from tau
rapprox  = sin( pi * tauhat/2)   %0.9888
%
% Kendal's Tau is Correlation between 
% rv sign(X1-X2) and sign(Y1-Y2)
signsx=[]; signsy=[];
for i = 1:n-1
    for j = i+1:n
      signsx = [signsx      sign(x(i) - x(j)) ];
      signsy = [signsy     sign(y(i) - y(j)) ];
      end
    end
corr(signsx', signsy')


