spear.m 2.4 KB

123456789101112131415161718192021222324252627282930313233343536373839404142434445464748495051525354555657585960616263646566676869707172737475767778798081828384858687888990919293949596979899100101102103104105106107108109110111112113114115116117118119120121122123124125126127128129130131132133
  1. %======================================================================
  2. % File: spear.m
  3. %
  4. %======================================================================
  5. function [r,t,p]=spear(x,y)
  6. %Syntax: [r,t,p]=spear(x,y)
  7. %__________________________
  8. %
  9. % Spearman's rank correlation coefficient.
  10. %
  11. % r is the Spearman's rank correlation coefficient.
  12. % t is the t-ratio of r.
  13. % p is the corresponding p-value.
  14. % x is the first data series (column).
  15. % y is the second data series, a matrix which may contain one or multiple
  16. % columns.
  17. %
  18. %
  19. % Reference:
  20. % Press W. H., Teukolsky S. A., Vetterling W. T., Flannery B. P.(1996):
  21. % Numerical Recipes in C, Cambridge University Press. Page 641.
  22. %
  23. %
  24. % Example:
  25. % x = [1 2 3 3 3]';
  26. % y = [1 2 2 4 3; rand(1,5)]';
  27. % [r,t,p] = spear(x,y)
  28. %
  29. %
  30. % Products Required:
  31. % Statistics Toolbox
  32. %
  33. % Alexandros Leontitsis
  34. % Department of Education
  35. % University of Ioannina
  36. % 45110- Dourouti
  37. % Ioannina
  38. % Greece
  39. %
  40. % University e-mail: me00743@cc.uoi.gr
  41. % Lifetime e-mail: leoaleq@yahoo.com
  42. % Homepage: http://www.geocities.com/CapeCanaveral/Lab/1421
  43. %
  44. % 3 Feb 2002.
  45. % x and y must have equal number of rows
  46. if size(x,1)~=size(y,1)
  47. error('x and y must have equal number of rows.');
  48. end
  49. % Find the data length
  50. N = length(x);
  51. % Get the ranks of x
  52. R = crank(x)';
  53. for i=1:size(y,2)
  54. % Get the ranks of y
  55. S = crank(y(:,i))';
  56. % Calculate the correlation coefficient
  57. r(i) = 1-6*sum((R-S).^2)/N/(N^2-1);
  58. end
  59. % Calculate the t statistic
  60. if r == 1 | r == -1
  61. t = r*inf;
  62. else
  63. t=r.*sqrt((N-2)./(1-r.^2));
  64. end
  65. % Calculate the p-values
  66. p=2*(1-tcdf(abs(t),N-2));
  67. function r=crank(x)
  68. %Syntax: r=crank(x)
  69. %__________________
  70. %
  71. % Assigns ranks on a data series x.
  72. %
  73. % r is the vector of the ranks
  74. % x is the data series. It must be sorted.
  75. %
  76. %
  77. % Reference:
  78. % Press W. H., Teukolsky S. A., Vetterling W. T., Flannery B. P.(1996):
  79. % Numerical Recipes in C, Cambridge University Press. Page 642.
  80. %
  81. %
  82. % Alexandros Leontitsis
  83. % Department of Education
  84. % University of Ioannina
  85. % 45110- Dourouti
  86. % Ioannina
  87. % Greece
  88. %
  89. % University e-mail: me00743@cc.uoi.gr
  90. % Lifetime e-mail: leoaleq@yahoo.com
  91. % Homepage: http://www.geocities.com/CapeCanaveral/Lab/1421
  92. %
  93. % 3 Feb 2002.
  94. u = unique(x);
  95. [xs,z1] = sort(x);
  96. [z1,z2] = sort(z1);
  97. r = (1:length(x))';
  98. r=r(z2);
  99. for i=1:length(u)
  100. s=find(u(i)==x);
  101. r(s,1) = mean(r(s));
  102. end