diff --git a/TheilSen.m b/TheilSen.m index ee1e3ac9b50a1b1a1b873ad88040063517191d80..13e5f8133f7fd0808941e06da7835f895e84af56 100644 --- a/TheilSen.m +++ b/TheilSen.m @@ -1,38 +1,40 @@ function [m, b] = TheilSen(data) -% Performs Theil-Sen robust linear regression on data +% Performs Theil-Sen robust linear regression on data. % -% [m b] = TheilSen(data) +% [m, b] = TheilSen(data) % % data: A MxD matrix with M observations. The first D-1 columns are the -% explanatory variables and the Dth column is the response such that -% data = [x1, x2, ..., x(D-1), y]; +% explanatory variables and the Dth column is the response such that +% data = [x1, x2, ..., x(D-1), y]; % % m: Estimated slope of each explanatory variable with respect to the -% response varuable. Therefore, m will be a vector of D-1 slopes. +% response varuable. Therefore, m will be a vector of D-1 slopes. % b: Estimated offsets. % % -% EXAMPLE: +% EXAMPLE % ------- -% n=100; -% outN=round(0.2*n); -% noise = randn(n,2)*0.1; noise(randi(n*2,[outN 1]))=randn(outN,1)*5; -% data = [linspace(0,10,n)' linspace(0,5,n)'] + noise; -% Bhat = [ones(n,1) data(:,1)]\data(:,2); +% n = 100; +% outN = round(0.2 * n); +% noise = randn(n, 2) * 0.1; +% noise(randi(n * 2, [outN, 1])) = randn(outN, 1) * 5; +% data = [linspace(0, 10, n)', linspace(0,5,n)'] + noise; +% Bhat = [ones(n, 1), data(:, 1)] \ data(:, 2); % [m, b] = TheilSen(data); -% plims = [min(data(:,1)) max(data(:,1))]'; -% figure -% plot(data(:,1),data(:,2),'k.',... -% plims,plims*m+b,'-r',... -% plims,plims*Bhat(2)+Bhat(1),'-b','linewidth',2) -% legend('Data','TheilSen','Least Squares','location','NW') -% title(sprintf('Acual Slope = %2.3g, LS est = %2.3g, TS est = %2.3g',[0.5 Bhat(2) m])) -% +% plims = [min(data(:, 1)), max(data(:, 1))]'; +% figure() +% plot(data(:, 1), data(:, 2), 'k.', ... +% plims, plims * m + b, '-r',... +% plims, plims * Bhat(2) + Bhat(1), '-b', 'linewidth', 2) +% legend('Data', 'Theil-Sen', 'Least Squares', 'location', 'NW') +% title(sprintf('Acual Slope = %2.3g, LS est. = %2.3g, TS est. = %2.3g', ... +% [0.5 Bhat(2) m])) +% % % Source: % Gilbert, Richard O. (1987), "6.5 Sen's Nonparametric Estimator of % Slope", Statistical Methods for Environmental Pollution Monitoring, -% John Wiley and Sons, pp. 217�219, ISBN 978-0-471-28878-7 +% John Wiley and Sons, pp. 217-219, ISBN 978-0-471-28878-7 % % % @@ -41,42 +43,44 @@ function [m, b] = TheilSen(data) % - updated help % - speed increase for 2D case % -% sz = size(data); -if length(sz)~=2 || sz(1)<2 +if length(sz) ~= 2 || sz(1) < 2 error('Expecting MxD data matrix with at least 2 observations.') end -if sz(2)==2 % normal 2-D case +if sz(2) == 2 % normal 2D case C = nan(sz(1)); - for i=1:sz(1) + + for i = 1:sz(1) % accumulate slopes - C(i,i:end) = (data(i,2)-data(i:end,2))./(data(i,1) - data(i:end,1)); + C(i, i:end) = (data(i, 2) - data(i:end, 2)) ./ ... + (data(i, 1) - data(i:end, 1)); end - m = nanmedian(C(:)); % calculate slope estimate + + m = nanmedian(C(:)); % calculate slope estimate - if nargout==2 - b = nanmedian(data(:,2)-m*data(:,1)); % calculate intercept if requested + if nargout == 2 + % calculate intercept if requested + b = nanmedian(data(:, 2) - m * data(:, 1)); end -else % other cases - C = nan(sz(1),sz(2)-1,sz(1)); - for i=1:sz(1) - C(:,:,i) = bsxfun( @rdivide,data(i,end)-data(:,end),... - bsxfun(@minus,data(i,1:end-1),data(:,1:end-1)) ); % accumulate slopes +else + C = nan(sz(1), sz(2)-1, sz(1)); + + for i = 1:sz(1) + % accumulate slopes + C(:, :, i) = bsxfun( @rdivide, data(i, end) - data(:, end), ... + bsxfun(@minus, data(i, 1:end-1), data(:, 1:end-1)) ); end - Cprm = reshape( permute(C,[1 3 2]), [],size(C,2),1 ); % stack layers of C to 2D - m = nanmedian(Cprm,1); % calculate slope estimate + + % stack layers of C to 2D + Cprm = reshape( permute(C, [1, 3, 2]), [], size(C, 2), 1 ); + m = nanmedian(Cprm, 1); % calculate slope estimate - if nargout==2 + if nargout == 2 % calculate all intercepts if requested - b = nanmedian( bsxfun(@minus,data(:,end),bsxfun(@times,m,data(:,1:end-1))) ); + b = nanmedian( bsxfun(@minus, data(:, end), ... + bsxfun(@times, m, data(:, 1:end-1))) ); end - end - - - - - diff --git a/license.txt b/license.txt index 62fe911adfcebafe36a5797aaf2f81004fe007df..bfb82553713528dd31951f19f541b780f2d6d848 100644 --- a/license.txt +++ b/license.txt @@ -5,11 +5,11 @@ Redistribution and use in source and binary forms, with or without modification, are permitted provided that the following conditions are met: - * Redistributions of source code must retain the above copyright - notice, this list of conditions and the following disclaimer. - * Redistributions in binary form must reproduce the above copyright - notice, this list of conditions and the following disclaimer in - the documentation and/or other materials provided with the distribution +* Redistributions of source code must retain the above copyright + notice, this list of conditions and the following disclaimer. +* Redistributions in binary form must reproduce the above copyright + notice, this list of conditions and the following disclaimer in + the documentation and/or other materials provided with the distribution THIS SOFTWARE IS PROVIDED BY THE COPYRIGHT HOLDERS AND CONTRIBUTORS "AS IS" AND ANY EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED TO, THE