rmtve.m

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

function [X1, T1X, C1X, error1]=RMTVE(X, max_fits, max_points)
% % RMTVE: Calculates the Reweighted Trimmed Minimum Vulume Ellipsoid (RMTVE) of a data set X.
% % 
% % Syntax;
% %
% % [X1,T1X,C1X]=RMTVE(X, max_fits, max_points);
% %
% % **********************************************************************
% % 
% % Description
% % 
% % Calculates the Reweighted Trimmed Minimum Vulume Ellipsoid (RMTVE)
% % of a data set X.
% %
% % This program is a modification of RMVE.  It has been modified to
% % trim the input data sets and trim the number of combinations of
% % line 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
% %
% % X is the matrix of the data set.
% %
% % 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
% %
% % X1 is the data contained in the RTMVE.
% %
% % T1X is the center of the RTMVE.
% %
% % C1X is the covariance matrix of the RTMVE.
% %
% % error1 is 1 if there is an error otherwise it is 0.
% %
% % **********************************************************************
% %
% % 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
% %
% % modified 14 February    2008    Trimmed the input data arrays.
% %                                 Updated comments.
% %                                 Improved the error handling and default
% %                                 values.
% % 
% % 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: MTVE, MVE, LMTSpolor, LMTSpol, LMSpol, LMTSreg, LMSreg, LMTSregor, LMSregor
% % 

if nargin < 1 || isempty(X) || ~isnumeric(X)
    warning('Not enough input arguments.  Return empty array.');
    error1=1;
    n=1;
    y=1;
    X1=[];
    T1X=[];
    C1X=[];
else

    pp=size(X,2);

    % If not input, set the maximum number of fits
    if nargin < 2  || 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 < 3 || 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
    
    
    % Calculate the Minimum Trimmed Vulume Ellipsoid (MTVE) 
    % of a data set X.
    [x50,x,TX,CX]=MTVE(X, max_fits, max_points);

    % n is the length of the data set, p is its dimension
    [n p]=size(X);

    % Chapter 7, eq. 1.28
    j=0;
    for i=1:n
        if (X(i,:)-TX)*inv(CX)*(X(i,:)-TX)'<=chi2inv(0.975,p)
            j=j+1;
            X1(j,:)=X(i,:);
        end
    end

    % Chapter 7, eq. 1.29
    T1X=mean(X1);

    % Chapter 7, eq. 1.30
    C1X=cov(X1);

end