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

Put example into its own script file.

parent 1df99aeb
No related branches found
No related tags found
No related merge requests found
......@@ -12,25 +12,6 @@ function [m, b] = TheilSen(data)
% 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', '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,
......
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', '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]))
0% Loading or .
You are about to add 0 people to the discussion. Proceed with caution.
Please register or to comment