function [r,t,p]=spear(x,y)
%Syntax: [r,t,p]=spear(x,y)
%__________________________
%
% Spearman's rank correalation coefficient.
%
% r is the Spearman's rank correlation coefficient.
% t is the t-ratio of r.
% p is the corresponding p-value.
% x is the first data series (column).
% y is the second data series, a matrix which may contain one or multiple
% columns.
%
%
% Reference:
% Press W. H., Teukolsky S. A., Vetterling W. T., Flannery B. P.(1996):
% Numerical Recipes in C, Cambridge University Press. Page 641.
%
%
% Example:
% x = [1 2 3 3 3]';
% y = [1 2 2 4 3; rand(1,5)]';
% [r,t,p] = spear(x,y)
%
%
% Products Required:
% Statistics Toolbox
%
% Alexandros Leontitsis
% Department of Education
% University of Ioannina
% 45110- Dourouti
% Ioannina
% Greece
% 
% University e-mail: me00743@cc.uoi.gr
% Lifetime e-mail: leoaleq@yahoo.com
% Homepage: http://www.geocities.com/CapeCanaveral/Lab/1421
% 
% 3 Feb 2002.

for i=1:size(y,2)

% quito NaN:
x=x(:);
yi=y(:,i);

H=find(x>=0 & yi>=0);
x=x(H);
yi=yi(H);

% x and y must have equal number of rows
if size(x,1)~=size(yi,1)
error('x and y must have equal number of rows.');
end

% Find the data length
N = length(x);

% Get the ranks of x
R = crank(x)';

% Get the ranks of y
S = crank(yi)';

% Calculate the correlation coefficient
r(i) = 1-6*sum((R-S).^2)/N/(N^2-1);

end

% Calculate the t statistic
if r == 1 | r == -1
t = r*inf;
else
t=r.*sqrt((N-2)./(1-r.^2));
end

% Calculate the p-values
p=2*(1-tcdf(abs(t),N-2));





