Improved documentation and example code.

.gitignore

+*.m~

README.md

-# Theil-Sen-Matlab
+# Theil-Sen estimator for Matlab
 
 ## Description
 
@@ -6,16 +6,18 @@ A fast, stand-alone Theil-Sen estimator for robust regression in Matlab.
 
 ### Theil-Sen estimator
 
-A [Theil-Sen estimator](https://en.wikipedia.org/wiki/Theil%E2%80%93Sen_estimator) provides robust linear regression in the 2D plane.
+A [Theil-Sen estimator](https://en.wikipedia.org/wiki/Theil%E2%80%93Sen_estimator) provides robust linear regression in the 2D plane:
 The resulting estimates of slope and intercept are relatively insensitive to outliers.
 
+FIXME: Mention main idea (data point pairs), and resulting complexity.
+
 ### No dependency on Matlab toolboxes
 
 This code is based on [Theil-Sen Robust Linear Regression](https://mathworks.com/matlabcentral/fileexchange/48294-theil-sen-robust-linear-regression), version 1.2.0.0, by [Zachary Danziger](https://mathworks.com/matlabcentral/profile/authors/1044524).
 It was modified mainly by re-defining its call to the `nanmedian()` function, which creates a dependency on the (commercially licensed) [Statistics Toolbox](https://mathworks.com/products/statistics.html).
 
 There are several other implementations of the Theil-Sen estimator on Mathworks's [File Exchange](https://mathworks.com/matlabcentral/fileexchange).
-Unfortunately all were found to depend on the Statistics Toolbox, [except one one](https://mathworks.com/matlabcentral/fileexchange/43135-regression-utilities), but that implementation was judged to be slower and less versatile.
+Unfortunately all were found to depend on the Statistics Toolbox, [except one](https://mathworks.com/matlabcentral/fileexchange/43135-regression-utilities), which was judged to be slower and less versatile.
 
 ## Installation
 
@@ -31,21 +33,27 @@ TODO: Add examples, and show some expected output.
     - [ ] Mention original use case.
     - [ ] Add usage example and plot.
 - [ ] Add a drop-in for `nanmean()` to become independent of Statistics Toolbox.
-- [ ] Adjust formatting to my taste.
-- [ ] Adapt license file.
+- [ ] Change input interface to match Matlab's `robustfit()`.
+- [ ] Add more checks for input sanity.
+
+
+## Changelog
+
+- October 2014 by Z. Danziger: Original version.
+- September 2015 by Z. Danziger: Updated help, speed increase for 2D case
+- March 2022 by J. Keyser: Adjusted formatting, added documentation,
 
 ## Contributing and project status
 
 This project is relatively unmaintained, and only shared as-is, in the hope to be helpful.
-The authors do not accept any liability.
-If you find a bug, feel free to the author(s) know (see below).
+If you find a bug, feel free to let the author(s) know.
 Feature requests should be directed to the original author.
 
 ## Authors
 
-1. Zachary Danziger, principal author ([Matlab profile](https://de.mathworks.com/matlabcentral/profile/authors/1044524), [Lab webpage](https://anil.fiu.edu/))
-2. Johannes Keyser [JLU GitLab profile](@jkeyser)
+1. Zachary Danziger, original author ([Matlab profile](https://de.mathworks.com/matlabcentral/profile/authors/1044524), [Lab webpage](https://anil.fiu.edu/))
+2. Johannes Keyser
 
 ## License
 
-Please see [license.txt](license.txt).
+[Simplified BSD License](https://en.wikipedia.org/wiki/BSD_licenses#2-clause_license_(%22Simplified_BSD_License%22_or_%22FreeBSD_License%22)), see [license.txt](license.txt).

TheilSen.m

@@ -3,27 +3,30 @@ function [m, b] = TheilSen(data)
 %
 % [m, b] = TheilSen(data)
 %
+% INPUT
 %   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];
 %
+% OUTPUT
 %   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
+%   See accompanying file example.m.
 %
-% Source:
-%   Gilbert, Richard O. (1987), "6.5 Sen's Nonparametric Estimator of
-%   Slope", Statistical Methods for Environmental Pollution Monitoring,
+% REFERENCE
+%   - 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
 %
+% AUTHORS
+%   2014-2015 Zachary Danziger
+%   2022 Johannes Keyser
 %
-%
-% %%% Z. Danziger October 2014 %%%
-% edits Z. Danziger September 2015:
-%	- updated help
-%   - speed increase for 2D case
-%
+% LICENSE
+%   Simplified BSD license, see accompanying file license.txt.
 
 sz = size(data);
 if length(sz) ~= 2 || sz(1) < 2

example.m

-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))]';
+% This example demonstrates the use of TheilSen.m with simulated data.
+clearvars *
+
+% Choose total sample size and fraction of outliers:
+N_total = 20;
+N_outlr = round(0.2 * N_total);
+
+% Simulate a linear data set
+true_b0 = -2;
+true_b1 = 10;
+data_x = linspace(0, 1, N_total)';
+data_y = true_b0 + true_b1 .* data_x;
+
+% add some "usual" Gaussian noise to both dimensions
+SDx_usual = 0.1;
+SDy_usual = 2.0;
+data_x = data_x + randn(size(data_x)) * SDx_usual;
+data_y = data_y + randn(size(data_y)) * SDy_usual;
+
+% simulate "outliers" as signal-dependent, positive noise in y dimension
+SDy_outlr = 20 * SDy_usual;
+outlr_idx = randperm(N_total, N_outlr);  % randomly select data points
+outlr_y = data_y(outlr_idx);  % start with original value
+outlr_y = abs(outlr_y) + 3 * SDy_usual;  % set to unlikely high value
+outlr_y = outlr_y + abs(randn(N_outlr, 1) * SDy_outlr);  % over-dispersion
+outlr_y = outlr_y .* data_x(outlr_idx);  % scale by x magnitude
+data_y(outlr_idx) = outlr_y;
+
+% estimate least squares parameters
+est_ls = [ones(N_total, 1), data_x] \ data_y;
+
+% estimate Theil-Sen parameters
+[m, b] = TheilSen([data_x, data_y]);
+est_ts = [b, m];
+
+% plot everything
 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]))
+plims = [min(data_x), max(data_x)]';
+normal_idx = setdiff(1:N_total, outlr_idx);
+plot(data_x(normal_idx), data_y(normal_idx), 'ko', ...
+     data_x(outlr_idx), data_y(outlr_idx), 'rx', ...
+     plims, plims * est_ls(2) + est_ls(1), '-r', ...
+     plims, plims * m + b, 'b-', ...
+     'linewidth', 2)
+legend('Normal data', 'Outlier', 'Least Squares', 'Theil-Sen', 'location', 'NW')
+title(sprintf('True: [%.3g, %.3g], LS: [%.3g, %.3g], TS: [%.3g, %.3g]', ...
+      [true_b0, true_b1, est_ls(1), est_ls(2), est_ts(1), est_ts(2)]))
\ No newline at end of file