12345678910111213141516171819202122232425262728293031323334353637383940414243444546474849505152535455565758596061626364656667686970717273 |
- function r = untiedrank(x)
- % function r = untiedrank(x)
- %
- % Similar to tiedrank, but arbitrarily breaks ties (with consistency each
- % time called)
- % reset random number generator to same start
- % RandStream.setDefaultStream(RandStream('mrg32k3a','Seed',10));
- RandStream.setGlobalStream(RandStream('mrg32k3a','Seed',10));
- if isvector(x)
- r = tr(x);
- else
- if isa(x,'single')
- outclass = 'single';
- else
- outclass = 'double';
- end
- % Operate on each column vector of the input (possibly > 2 dimensional)
- sz = size(x);
- ncols = sz(2:end); % for 2x3x4, ncols will be [3 4]
- r = zeros(sz,outclass);
- for j=1:prod(ncols)
- r(:,j)= tr(x(:,j));
- end
- end
- % --------------------------------
- function r = tr(x)
- %TR Local untiedrank function to compute results for one column
- % Sort, then leave the NaNs (which are sorted to the end) alone
- [sx, rowidx] = sort(x(:));
- numNaNs = sum(isnan(x));
- xLen = numel(x) - numNaNs;
- % Use ranks counting from low end
- ranks = [1:xLen NaN(1,numNaNs)]';
- if isa(x,'single')
- ranks = single(ranks);
- end
- % "randomly" break ties. Avoid using diff(sx) here in case there are infs.
- ties = (sx(1:xLen-1) == sx(2:xLen));
- tieloc = [find(ties); xLen+2];
- maxTies = numel(tieloc);
- tiecount = 1;
- while (tiecount < maxTies)
- tiestart = tieloc(tiecount);
- ntied = 2;
- while(tieloc(tiecount+1) == tieloc(tiecount)+1)
- tiecount = tiecount+1;
- ntied = ntied+1;
- end
-
- % Compute mean of tied ranks
- % ranks(tiestart:tiestart+ntied-1) = sum(ranks(tiestart:tiestart+ntied-1)) / ntied;
- % "randomly" reassign ties
- temp = ranks(tiestart:tiestart+ntied-1);
- ranks(tiestart:tiestart+ntied-1) = temp(randperm(ntied));
-
- tiecount = tiecount + 1;
- end
- % Broadcast the ranks back out, including NaN where required.
- r(rowidx) = ranks;
- r = reshape(r,size(x));
|