--- scripts/statistics/base/ranks.m 2015-12-20 12:45:06.832851328 +1100 +++ /ranks_20160108.m 2016-01-08 16:46:11.718096775 +1100 @@ -1,4 +1,5 @@ ## Copyright (C) 1995-2015 Kurt Hornik +## Copyright (C) 2012 Dave Goel ## ## This file is part of Octave. ## @@ -17,24 +18,33 @@ ## . ## -*- texinfo -*- -## @deftypefn {} {} ranks (@var{x}, @var{dim}) -## Return the ranks of @var{x} along the first non-singleton dimension -## adjusted for ties. +## @deftypefn {} {} ranks (@var{x}) +## @deftypefnx {} {} ranks (@var{x}, @var{dim}) +## @deftypefnx {} {} ranks (@var{x}, @var{dim}, @var{rtype}) ## -## If the optional argument @var{dim} is given, operate along this dimension. +## Return the ranks (in the sense of order statistics) of @var{x} along +## the first non-singleton dimension +## adjusted for ties. If the argument @var{dim} is given and non-empty, +## operate along this dimension instead. +## +## The parameter @var{rtype} determines the way ties are handled: +## 0 or "fractional" (the default) for fractional, which gives ranks +## like (1, 2.5, 2.5, 4); +## 1 or "competition" for competition ranking (Example: [1, 2, 2, 4], +## if @var{x} is [11 12 12 14]; +## 2 or "modified" for modified competition ranking (1, 3, 3, 4); +## 3 or "ordinal" for ordinal ranking (1, 2, 3, 4); +## 4 or "revordinal" for reverse ordinal (1, 3, 2, 4)), and +## 5 or "dense" for dense ranking (1, 2, 2, 3). ## @seealso{spearman, kendall} ## @end deftypefn -## Author: KH -## Description: Compute ranks - -## This code was rather ugly, since it didn't use sort due to the -## fact of how to deal with ties. Now it does use sort and its -## even uglier!!! At least it handles NDArrays.. +## Authors: Dave Goel , KH , +## Description: Compute ranks using a choice of methods when resolving ties. -function y = ranks (x, dim) +function y = ranks (x, dim, rtype) - if (nargin != 1 && nargin != 2) + if (nargin < 1) || (nargin > 3); print_usage (); endif @@ -42,9 +52,14 @@ error ("ranks: X must be a numeric vector or matrix"); endif + if (nargin < 3) || ( isempty (rtype)); + rtype = 0; + endif + nd = ndims (x); sz = size (x); - if (nargin != 2) + + if (nargin < 2) || isempty (dim); ## Find the first non-singleton dimension. (dim = find (sz > 1, 1)) || (dim = 1); else @@ -55,35 +70,87 @@ endif if (sz(dim) == 1) - y = ones (sz); + y = ones (sz); ## Dimension dim is singleton, so all are ranked first. else - ## The algorithm works only on dim = 1, so permute if necesary. + ## The algorithm works only on dim = 1, so permute if necessary. if (dim != 1) perm = [1 : nd]; perm(1) = dim; perm(dim) = 1; x = permute (x, perm); + sz = size (x); endif - sz = size (x); - infvec = -Inf ([1, sz(2 : end)]); - [xs, xi] = sort (x); - eq_el = find (diff ([xs; infvec]) == 0); - if (isempty (eq_el)) - [eq_el, y] = sort (xi); - else - runs = setdiff (eq_el, eq_el+1); - len = diff (find (diff ([Inf; eq_el; -Inf]) != 1)) + 1; - [eq_el, y] = sort (xi); - for i = 1 : length (runs) - y (xi (runs (i) + [0:(len(i)-1)]) + floor (runs (i) ./ sz(1)) - * sz(1)) = eq_el(runs(i)) + (len(i) - 1) / 2; - endfor - endif + + switch rtype; + case {4, "revordinal"}; + [sx ids] = sort (x, 'descend'); + ids = flipdim (ids, 1); + otherwise + [sx ids] = sort (x); ## sx is sorted x. + endswitch + lin = cumsum (ones (size (x)), 1); ## A linearly increasing array. + + switch rtype; + case {0, "fractional"}; + lin = (_competition (lin, sx, sz)+_modified (lin, sx, sz)) / 2; + case {1, "competition"}; + lin = _competition (lin, sx, sz); + case {2, "modified"}; + lin = _modified (lin, sx, sz); + case {3, "ordinal"}; + ## no processing needed here. + case {4, "revordinal"}; + ## no processing needed here. + case {5, "dense"}; + lin = _dense (lin, sx, sz); + otherwise + if (!ischar (rtype)) + rtype = num2str (rtype); + end + error ("ranks: Illegal RTYPE '%s'", rtype); + endswitch + y = NaN (size (lin)); + + + ## Offsets to map indices into each column to indices into the linear array. + idf = zeros (sz); + idf(1, :) = 0:sz(1): (numel (ids)-1); + idf(:, :) = repmat (idf (1, :), [sz(1), ones (1, length (sz) - 1)]); + y(ids + idf) = lin; + if (dim != 1) y = permute (y, perm); endif endif +endfunction + + + +function linnew = _dense (lin, sx, sz) + infvec = -Inf * ones ([1, sz(2 : end)]); + fnewp = logical (diff ([infvec; sx])); + linnew = cumsum (fnewp, 1); +endfunction + +function linnew = _competition (lin, sx, sz) + ## Stop increasing lin when sx does not increase. Same as before + ## otherwise. + infvec = -Inf * ones ([1, sz(2 : end)]); + fnewp = find (diff ([infvec; sx])); + linnew = zeros (size (lin)); + linnew(fnewp) = lin(fnewp); + linnew = cummax (linnew, 1); +endfunction + +function linnew = _modified (lin, sx, sz) + ## Traverse lin backwards. Stop decreasing it when sx doesn't + ## decrease. + infvec = Inf * ones ([1, sz(2 : end)]); + fnewp = find (diff ([sx; infvec])); + linnew = Inf (size (lin)); + linnew(fnewp) = lin(fnewp); + linnew = flipdim (cummin (flipdim (linnew, 1)), 1); endfunction @@ -96,7 +163,7 @@ ## Test input validation %!error ranks () -%!error ranks (1, 2, 3) +%!error ranks (1, 2, 3, 4) %!error ranks ({1, 2}) %!error ranks (['A'; 'B']) %!error ranks (1, 1.5)