Skip to content
Snippets Groups Projects
Commit 34a11524 authored by Keyser, Johannes's avatar Keyser, Johannes
Browse files

Make internal notation of terms explicit/simpler

parent 2e8626e4
No related branches found
No related tags found
No related merge requests found
function [m, b] = TheilSen(data) function [b1, b0] = TheilSen(data)
% Performs Theil-Sen robust linear regression on data. % Performs Theil-Sen robust linear regression on data.
% %
% [m, b] = TheilSen(data) % [b1, b0] = TheilSen(data)
% %
% INPUT % INPUT
% data: A MxD matrix with M observations. The first D-1 columns are the % data: A Num_Obs x Num_Dim matrix with Num_Obs observations.
% explanatory variables and the Dth column is the response such that % The first Num_Dim - 1 columns are the explanatory variables and the
% data = [x1, x2, ..., x(D-1), y]; % last column is the response such that
% data = [x1, x2, ..., x(Num_Dim - 1), y];
% %
% OUTPUT % OUTPUT
% m: Estimated slope of each explanatory variable with respect to the % b1: Estimated slope of each explanatory variable with respect to the
% response varuable. Therefore, m will be a vector of D-1 slopes. % response variable. Therefore, b1 will be a vector of Num_Dim - 1 slopes.
% b: Estimated offsets. % b0: Estimated offsets.
% %
% EXAMPLE % EXAMPLE
% See accompanying file example.m. % See accompanying file example.m.
...@@ -22,37 +23,45 @@ function [m, b] = TheilSen(data) ...@@ -22,37 +23,45 @@ function [m, b] = TheilSen(data)
% 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
% %
% AUTHORS % AUTHORS
% 2014-2015 Zachary Danziger % 2014-15 Zachary Danziger
% 2022 Johannes Keyser % 2022 Johannes Keyser
% %
% LICENSE % LICENSE
% Simplified BSD license, see accompanying file license.txt. % BSD 2-clause "simplified" license, see accompanying file license.txt.
sz = size(data); sz = size(data);
if length(sz) ~= 2 || sz(1) < 2
error('Expecting MxD data matrix with at least 2 observations.') if length(sz) ~= 2
error('Expecting a 2D data matrix Num_Obs x Num_Dim.')
end
Num_Obs = sz(1); % number of observations
Num_Dim = sz(2); % number of dimensions
if Num_Obs < 2
error('Expecting a data matrix Obs x Dim with at least 2 observations.')
end end
if sz(2) == 2 % normal 2D case if Num_Dim == 2 % normal 2D case
C = nan(sz(1)); C = nan(Num_Obs);
for i = 1:sz(1) for i = 1:Num_Obs
% accumulate slopes % accumulate slopes
C(i, i:end) = (data(i, 2) - data(i:end, 2)) ./ ... C(i, i:end) = (data(i, 2) - data(i:end, 2)) ./ ...
(data(i, 1) - data(i:end, 1)); (data(i, 1) - data(i:end, 1));
end end
m = nanmedian(C(:)); % calculate slope estimate b1 = nanmedian(C(:)); % calculate slope estimate
if nargout == 2 if nargout == 2
% calculate intercept if requested % calculate intercept if requested
b = nanmedian(data(:, 2) - m * data(:, 1)); b0 = nanmedian(data(:, 2) - b1 * data(:, 1));
end end
else else
C = nan(sz(1), sz(2)-1, sz(1)); C = nan(Num_Obs, Num_Dim-1, Num_Obs);
for i = 1:sz(1) for i = 1:Num_Obs
% accumulate slopes % accumulate slopes
C(:, :, i) = bsxfun( @rdivide, data(i, end) - data(:, end), ... C(:, :, i) = bsxfun( @rdivide, data(i, end) - data(:, end), ...
bsxfun(@minus, data(i, 1:end-1), data(:, 1:end-1)) ); bsxfun(@minus, data(i, 1:end-1), data(:, 1:end-1)) );
...@@ -60,11 +69,11 @@ else ...@@ -60,11 +69,11 @@ else
% stack layers of C to 2D % stack layers of C to 2D
Cprm = reshape( permute(C, [1, 3, 2]), [], size(C, 2), 1 ); Cprm = reshape( permute(C, [1, 3, 2]), [], size(C, 2), 1 );
m = nanmedian(Cprm, 1); % calculate slope estimate b1 = nanmedian(Cprm, 1); % calculate slope estimate
if nargout == 2 if nargout == 2
% calculate all intercepts if requested % calculate all intercepts if requested
b = nanmedian( bsxfun(@minus, data(:, end), ... b0 = nanmedian( bsxfun(@minus, data(:, end), ...
bsxfun(@times, m, data(:, 1:end-1))) ); bsxfun(@times, b1, data(:, 1:end-1))) );
end end
end end
\ No newline at end of file
0% Loading or .
You are about to add 0 people to the discussion. Proceed with caution.
Please register or to comment