lmtspol.m

This toolbox implements the same methods on small dadta sets and imlements a trimming method using a

function [LMSout, blms, Rsq, error1]=LMTSpol(y, p, x, max_fits, max_points)
% % LMTSpol: Estimates the best fit polynomial using a random trimming and least median squares
% % 
% % Syntax;
% % 
% % [LMSout, blms, Rsq]=LMTSpol(y, p, x, max_fits, max_points);
% % 
% % **********************************************************************
% % 
% % Description
% % 
% % Calculates the Least Median Trimmed Squares (LMS) polynomial 
% % regression parameters and output. 
% % 
% % This program is a modification of LMSpol.  It has been modified to 
% % randomly trim the input data sets and trim the number of combinations 
% % of polynomial fits that are processed.  The trimming allows the program
% % to accomodate large data sets.
% % 
% % This program performs the Least Median Trimmed Squares Robust
% % Regression for simple or multiple columns of data and outputs the
% % regression parameters. 
% % 
% % **********************************************************************
% % 
% % Input Variable Description
% % 
% % y is the vector of the dependent variable.
% % 
% % p is the order of the polynomial.
% % 
% % x is the vector of the independent variable.
% % 
% % max_fits is the number of best fit pairs of data. 
% %     The maximum value is 10000.
% %     The default value is 1000 or the largest value allowed.  
% % 
% % max_points is the number of data points for curve fitting.  
% %     The maximum value is 100000.  
% %     The default value is 100000 or the largest value allowed.
% % 
% %
% % **********************************************************************
% % 
% % Output Variable Description
% %
% % LMSout is the LMS estimated values vector.
% % 
% % blms is the LMS [intercept slopes] vector of the form
% %       y=p0+p1*x+p2*x^2+p3*x3 ... pn*x^n;
% %       blms=[p0 p1 p2 p3 p4 p5 ... pn];
% % 
% % Rsq is the R-squared regression coefficient error estimate of the fit.
% % 
% % error1 is 1 if there is an error otherwise it is 0.
% %
% % **********************************************************************
% 
% Example='1';
% % Establish an exact solution (xe, ye)
% 
% p=3; % The polynomial will be third order.
% xe=(1:100)';
% ye=1+3.*xe+2.*xe.^2+3.*xe.^3;
%  
% % Create a noisy test signal
% x=xe+randn(size(xe));
% y=1+randn(size(xe))+3.*(xe+randn(size(xe)))+2.*(xe+randn(size(xe))).^2+3.*(xe+randn(size(xe))).^3;
% 
% 
% 
% % Perform the robust median trimmed squares linear regression
% max_fits=100;
% max_points=500;
%
% % Outlier data points form a line with opposite slope
% % randomly select pcnt of the data points to be outliers
% pcnt=40;
% [ndraw]=rand_int(1, length(xe), pcnt/100*length(xe), 1, 1);
% x(ndraw)=ndraw+randn(size(ndraw));
% y(ndraw)=-1-randn(size(ndraw))-3.*(x(ndraw)+randn(size(ndraw)))-2.*(x(ndraw)+randn(size(ndraw))).^2-3.*(x(ndraw)+randn(size(ndraw))).^3;
% 
% % Sort data in ascending order
% [x, ix]=sort(x);
% y=y(ix);
%  
% % Perform the Robust Fit
% [LMSout,blms,Rsq]=LMTSpol(y, p, x, 1000, 5000);
% yest=polyval(flipud(blms), x);
% 
% % Perform the regular fit
% pcoefs = polyfit(x, y, p);
% yest2=polyval(pcoefs, x);
% 
% % Plot the Results
%
% figure(1); 
% plot(x, y', 'linestyle', 'none', 'marker', '.', 'markersize', 5, 'markeredgecolor', 'k');
% hold on; 
% plot(xe, ye', 'g');
% plot(x, yest', '--r');
% plot(x, yest2', ':b');
% legend({'Scattered Data', 'Exact Solution', 'Robust Solution', 'Regular Regression'});
% xlim([1 100]);
% title({[num2str(100-pcnt), '% of the data are good'], [num2str(pcnt), '% of the data are outliers']}, 'fontsize', 20);
% xlabel('x-axis', 'fontsize', 18);
% ylabel('y-axis', 'fontsize', 18);
% set(gca, 'fontsize', 14);
% 
% 
% 
% 
% % **********************************************************************
% % 
% % Reference:
% % Rousseeuw PJ, Leroy AM (1987): Robust regression and outlier detection.
% % Wiley.
% % 
% % **********************************************************************
% % 
% % This program is originally the work of
% % 
% % Alexandros Leontitsis
% % Institute of Mathematics and Statistics
% % University of Kent at Canterbury
% % Canterbury
% % Kent, CT2 7NF
% % U.K.
%
% % University e-mail: [e-mail] (until December 2002)
% % Lifetime e-mail: [e-mail]
% % Homepage: http://www.geocities.com/CapeCanaveral/Lab/1421
% % 
% % Sep 3, 2001.
% % 
% % 
% %
% % **********************************************************************
% %
% %
% % This program was modified by Edward L. Zechmann
% %
% %     date  1 February    2008    Updated comments.
% %                                 Added LMTSreg.
% %                                 Updated the logic syntax.
% %                             
% % modified 14 February    2008    Updated comments.
% %
% % modified  2 December    2008    Fixed a bug in trimming the input data 
% %                                 arrays.  This fix improves accuracy 
% %                                 for data sets with less than 1000 
% %                                 points.  
% %
% %
% % 
% % **********************************************************************
% %
% % Feel free to modify this code.
% % 
% % See also: LMSpol, LMTSpolor, LMSpolor, LMTSreg, LMSreg, LMTSregor, LMSregor
% % 



if logical(nargin < 1) || isempty(y) || ~isnumeric(y)
   error('Not enough input arguments.');
else
   % y must be a column vector
   y=y(:);
   % n is the length of the data set
   n=length(y);
end

if logical(nargin < 2) || isempty(p) || ~isnumeric(p)
   % If p is omitted give it the value of 1
   p=1;
else
   % p must be a scalar
   if logical( numel(p) > 1 )
      error('p must be a scalar.');
   end   
   % p must be a non-negative integrer
   if logical(round(p)-p~=0) || logical(p < 0)
      error('p must be a non-negative integer');
   end
end

if logical(nargin < 3) || isempty(x) || ~isnumeric(x)
   % If x is omitted give it the values 1:n
   x=(1:n)';
else
   % x must be a column vector
   x=x(:);
   % x and y must have the same length
   if n~=size(x,1)
      error('x and y must have the same length.');
   end
end

if n <= p
   error('The polynomial order is too large for the data set.');
end

% Prepare the matrix X for regression.
X=zeros(n, p);
for ix=1:p;
   X(:,ix)=x.^ix;
end


pp=size(X,2);

% If not input, set the maximum number of fits
if nargin < 4 || isempty(max_fits) || ~isnumeric(max_fits)
    % default value of max_fits is 1000
    max_fits=min([1000, nchoosek(n, pp+1)]);
end

% make sure that max_fits does not exceed 10000
max_fits=min( [max_fits, nchoosek(n, pp+1), 10000]);

% If max_points is not an input, set the maximum number of points
% for the input arrays X and y to a reasonable value.
if nargin < 5 || isempty(max_points) || ~isnumeric(max_points)
    max_points=max([min([n, 100000]), max_fits*(pp+1)]);
end

if max_points < max_fits
    max_points=max_fits;
end


% Perform the LMTS regression
[LMSout, blms, Rsq, error1]=LMTSreg(y, X, max_fits, max_points);