untiedrank.m 1.8 KB

12345678910111213141516171819202122232425262728293031323334353637383940414243444546474849505152535455565758596061626364656667686970717273
  1. function r = untiedrank(x)
  2. % function r = untiedrank(x)
  3. %
  4. % Similar to tiedrank, but arbitrarily breaks ties (with consistency each
  5. % time called)
  6. % reset random number generator to same start
  7. % RandStream.setDefaultStream(RandStream('mrg32k3a','Seed',10));
  8. RandStream.setGlobalStream(RandStream('mrg32k3a','Seed',10));
  9. if isvector(x)
  10. r = tr(x);
  11. else
  12. if isa(x,'single')
  13. outclass = 'single';
  14. else
  15. outclass = 'double';
  16. end
  17. % Operate on each column vector of the input (possibly > 2 dimensional)
  18. sz = size(x);
  19. ncols = sz(2:end); % for 2x3x4, ncols will be [3 4]
  20. r = zeros(sz,outclass);
  21. for j=1:prod(ncols)
  22. r(:,j)= tr(x(:,j));
  23. end
  24. end
  25. % --------------------------------
  26. function r = tr(x)
  27. %TR Local untiedrank function to compute results for one column
  28. % Sort, then leave the NaNs (which are sorted to the end) alone
  29. [sx, rowidx] = sort(x(:));
  30. numNaNs = sum(isnan(x));
  31. xLen = numel(x) - numNaNs;
  32. % Use ranks counting from low end
  33. ranks = [1:xLen NaN(1,numNaNs)]';
  34. if isa(x,'single')
  35. ranks = single(ranks);
  36. end
  37. % "randomly" break ties. Avoid using diff(sx) here in case there are infs.
  38. ties = (sx(1:xLen-1) == sx(2:xLen));
  39. tieloc = [find(ties); xLen+2];
  40. maxTies = numel(tieloc);
  41. tiecount = 1;
  42. while (tiecount < maxTies)
  43. tiestart = tieloc(tiecount);
  44. ntied = 2;
  45. while(tieloc(tiecount+1) == tieloc(tiecount)+1)
  46. tiecount = tiecount+1;
  47. ntied = ntied+1;
  48. end
  49. % Compute mean of tied ranks
  50. % ranks(tiestart:tiestart+ntied-1) = sum(ranks(tiestart:tiestart+ntied-1)) / ntied;
  51. % "randomly" reassign ties
  52. temp = ranks(tiestart:tiestart+ntied-1);
  53. ranks(tiestart:tiestart+ntied-1) = temp(randperm(ntied));
  54. tiecount = tiecount + 1;
  55. end
  56. % Broadcast the ranks back out, including NaN where required.
  57. r(rowidx) = ranks;
  58. r = reshape(r,size(x));