% SPCORR calculate Spearman Rank-Order Correlation Coefficient
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% Description:
% 
% SPCORR calculates Spearman Rank-Order Correlation Coefficient
% 
% CALL:
%       [rs,pValue,TestRes] = spcorr(y1,y2,alpha)
% 
% INPUT: 
%       - y1 and y2 are data series
%       - alpha is the significance level. Default 0.05
% 
% OUTPUT:
%       - rs: Spearman rank-order correlation coefficient
%       - pValue: Probability associated with the calculated rs in the
%       hypothesis that the correlation between y1 and y2 is zero
%       - TestRes: True or false on the basis of the test results. The null
%       hypothesis for the test is that the two series y1 and y2 are
%       uncorrelated.
%         TestRes = 0 => Do not reject the null hypothesis at significance
%         level alpha. (pValue >= alpha)
%         TestRes = 1 => Reject the null hypothesis at significance level
%         alpha. (pValue < alpha)
% 
% NOTE: 
%       The statistic of Spearman rank-order correlation coefficient is
%       well approximated by a Student t distribution. Hypothesis test is
%       then based on such statistic.
% 
%   References:
%      [1] W. H. Press, S. A. Teukolsky, W. T. Vetterling, B. P. Flannery,
%      Numerical Recipes 3rd Edition: The Art of Scientific Computing,
%      Cambridge University Press; 3 edition (September 10, 2007)
% 
%
% L Ferraioli 06-12-2010
%
% $Id: spcorr.m,v 1.2 2011/07/06 14:42:26 luigi Exp $
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
function [rs,pValue,TestRes] = spcorr(y1,y2,alpha)
  
  % check size
  if size(y1)~=size(y2)
    error('y1 and y2 must have the same size')
  end
  
  % check input
  if isempty(alpha)
    alpha = 0.05;
  elseif alpha > 1
    alpha = alpha/100;
  end
    
  
  % calculate rank for y1 and y2
  [r1,s1] = utils.math.crank(y1);
  [r2,s2] = utils.math.crank(y2);
  
  n = numel(y1);
  
  dv = (r1-r2).^2;
  
  dd = sum(dv);
  
  en = n*n*n-n;
  
  sq1 = sqrt(1-s1/en);
  sq2 = sqrt(1-s2/en);
  
  % calculate Spearman rank-order correlation coefficient
  rs = (1 - 6*(dd + s1/12 + s2/12)/en)/(sq1*sq2);
  
  % transform rs in t which is according Student's distribution with n-2
  % degrees of freedom
  t = rs*sqrt((n-2)/(1-rs^2));
  
  % Indicated with f(x) the prob. distribution function (Student's t in the
  % present caase). The probability Pr(k <= t) is proportional to the
  % integral Int_(-t,t)[f(x)dx]. For a Student's distribution such integral
  % is represeted by the Beta incomplete function
  % Pr(k <= t) = betainc(df/(df+t^2),df/2,1/2,'upper')
  % As a consequence Pr(k > t) = 1 - Pr(k <= t) which in Matlab can be
  % effectively calculated as
  % Pr(k > t) = betainc(df/(df+t^2),df/2,1/2,'lower')
  % Which provides the probability of finding a value grater than t if our
  % variable k is distributed according to a Student's distribution with df
  % degrees of freedom. Such a value represent the required pValue for the
  % test in the hypothesis that the two input data series are uncorrelated
  fac = (rs+1)*(1-rs);
  if fac>0
    df = n-2;
    pValue = betainc(df/(df+t*t),0.5*df,0.5,'lower');
  else
    pValue = 0;
  end
  
  TestRes = (pValue < alpha);
  
  
end