From 913d03c637dddff040163c8dad74a0fb485514d1 Mon Sep 17 00:00:00 2001 From: Johannes Keyser <johannes.keyser@sport.uni-giessen.de> Date: Mon, 21 Mar 2022 19:13:08 +0100 Subject: [PATCH] add unchanged code+license from File Exchange --- TheilSen.m | 82 +++++++++++++++++++++++++++++++++++++++++++++++++++++ license.txt | 24 ++++++++++++++++ 2 files changed, 106 insertions(+) create mode 100644 TheilSen.m create mode 100644 license.txt diff --git a/TheilSen.m b/TheilSen.m new file mode 100644 index 0000000..ee1e3ac --- /dev/null +++ b/TheilSen.m @@ -0,0 +1,82 @@ +function [m, b] = TheilSen(data) +% Performs Theil-Sen robust linear regression on 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]; +% +% m: Estimated slope of each explanatory variable with respect to the +% response varuable. Therefore, m will be a vector of D-1 slopes. +% b: Estimated offsets. +% +% +% 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); +% [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])) +% +% +% 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 +% +% +% +% %%% Z. Danziger October 2014 %%% +% edits Z. Danziger September 2015: +% - updated help +% - speed increase for 2D case +% +% + +sz = size(data); +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 + C = nan(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)); + end + m = nanmedian(C(:)); % calculate slope estimate + + if nargout==2 + b = nanmedian(data(:,2)-m*data(:,1)); % calculate intercept if requested + 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 + 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 + + if nargout==2 + % calculate all intercepts if requested + b = nanmedian( bsxfun(@minus,data(:,end),bsxfun(@times,m,data(:,1:end-1))) ); + end + +end + + + + + diff --git a/license.txt b/license.txt new file mode 100644 index 0000000..62fe911 --- /dev/null +++ b/license.txt @@ -0,0 +1,24 @@ +Copyright (c) 2015, Zachary Danziger +All rights reserved. + +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 + +THIS SOFTWARE IS PROVIDED BY THE COPYRIGHT HOLDERS AND CONTRIBUTORS "AS IS" +AND ANY EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED TO, THE +IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR PURPOSE +ARE DISCLAIMED. IN NO EVENT SHALL THE COPYRIGHT OWNER OR CONTRIBUTORS BE +LIABLE FOR ANY DIRECT, INDIRECT, INCIDENTAL, SPECIAL, EXEMPLARY, OR +CONSEQUENTIAL DAMAGES (INCLUDING, BUT NOT LIMITED TO, PROCUREMENT OF +SUBSTITUTE GOODS OR SERVICES; LOSS OF USE, DATA, OR PROFITS; OR BUSINESS +INTERRUPTION) HOWEVER CAUSED AND ON ANY THEORY OF LIABILITY, WHETHER IN +CONTRACT, STRICT LIABILITY, OR TORT (INCLUDING NEGLIGENCE OR OTHERWISE) +ARISING IN ANY WAY OUT OF THE USE OF THIS SOFTWARE, EVEN IF ADVISED OF THE +POSSIBILITY OF SUCH DAMAGE. -- GitLab