mtve.m

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

function [x50,x,TX,CX]=MTVE(X, max_fits, max_points)
% % MTVE: Calculates the Minimum Trimmed Vulume Ellipsoid (MTVE) of a data set X.
% % 
% % Syntax;
% %
% % [x50,x,TX,CX]=MTVE(X, max_fits, max_points);
% %
% % **********************************************************************
% % 
% % Description
% % 
% % Calculates the Minimum Trimmed Vulume Ellipsoid (MTVE) of a data set X.
% %
% % This program is a modification of MVE.  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
% %
% % x50 is the 50% of the points that form the MTVE.
% %
% % x is the core subsample for the MTVE.
% %
% % TX is the center of the MTVE.
% %
% % CX is the inflated covariance matrix of the MTVE.
% %
% % **********************************************************************
% %
% % 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: MVE, LMSpolor, LMTSpol, LMSpol, LMTSreg, LMSreg, LMTSregor, LMSregor
% % 

flag=0;

if nargin < 1 || isempty(X) || ~isnumeric(X)
    warning('Not enough input arguments.  Return empty array.');
    flag=1;
else
    % X must be 2-dimensional
    if ndims(X) > 2
        warning('Invalid data set.  Return empty array.');
        flag=1;
    end

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

if n < p
    warning('You must give a larger data set.  Return empty array.');
    flag=1;
end

if flag == 1
    x50=[];
    x=[];
    TX=[];
    CX=[];
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


    % Program Modified Here
    % input data is trimmed
    % best fit combinations are trimmed
    y=X;
    [C, y, X, n, p]=LMS_trim(y, X, max_fits, max_points, 1);
    

    volmin=Inf;

    for i=1:size(C,1);
        
        for j=1:p+1
            A(j,:)=X(C(i,j),:);
        end
        
        if rank(A)==p

            %Chapter 7, eq. 1.23
            Cj=cov(A);

            %Chapter 7, eq. 1.24
            for j=1:n
                fact(j)=(X(j,:)-mean(A))*inv(Cj)*(X(j,:)-mean(A))';
            end
            mj2=median(fact);

            mj=sqrt(mj2);

            % Chapter 7, eq. 1.25 - The objective function
            vol=sqrt(det(Cj))*mj^(p-1);

            if vol<volmin
                volmin=vol;

                % Chapter 7, eq. 1.26 - MVE
                TX=mean(A);
                CX=mj2*Cj/chi2inv(0.5,p);

                % The core of the MVE
                x=A;

                % The 50% of the points contained in the MVE
                k=0;
                x50=[];
                for j=1:n
                    if fact(j)<=mj2
                        k=k+1;
                        x50(k,:)=X(j,:);
                    end
                end
            end
        end
    end

end