ECG-Kit 1.0

function result=mloclogist(x,varargin)

%MLOCLOGIST is an M-estimator of location with psi-function equal to
% (exp(x)-1)/(exp(x)+1) and with an auxiliary scale estimate.
% It is iteratively computed. It can resist 50% outliers. 
% If x is a matrix, the location estimate is computed on the columns of x. The
% result is then a row vector. If x is a row or a column vector,
% the output is a scalar.
%
% The estimator is introduced in:
%   Rousseeuw, P.J. and Verboven, S. (2002),
%   "Robust estimation in very small samples", 
%   Computational Statistics and Data Analysis, 40, 741-758.
%
% Required input argument:
%    x: either a data matrix with n observations in rows, p variables in columns
%       or a vector of length n.
%
% Optional input arguments:
%    k: number of iteration steps (default value = 50)
%  loc: a starting value for the location estimate
%       default value = 'median'
%       other possibilities : 'hl'/'mlochuber'/...
%  sca: an auxiliary scale estimate
%       default value = 'madc'
%       other possibilities: 'qn'/'adm'/'mscalelogist'/...
%
% I/O: result=mloclogist(x,'k',50,'loc','median','sca','madc')    
%
% Examples: result=mloclogist(x,'loc','mlochuber','sca','qn');
%           result=mloclogist(x,'sca','mscalelogist','k',10);
%
% This function is part of LIBRA: the Matlab Library for Robust Analysis,
% available at:  
%              http://wis.kuleuven.be/stat/robust.html
%
% Written by S. Verboven
%Revisions by N. Smets
% Last update 28/08/03 

[n,p]=size(x);
%
% initialization with defaults
%
counter=1;
sca='madc';
loc='median';
default=struct('k',50,'sca',sca,'loc',loc);
list=fieldnames(default);
options=default;
IN=length(list);
i=1;
%
%
if nargin>1
    %
    % placing inputfields in array of strings
    %
    for j=1:nargin-1  
        if rem(j,2)~=0
            chklist{i}= varargin{j};    
            i=i+1;
        end
    end 
    %
    % Checking which default parameters have to be changed
    % and keep them in the structure 'options'.
    %
    while counter<=IN 
        index=strmatch(list(counter,:),chklist,'exact');
        if ~isempty(index) % in case of similarity
            for j=1:nargin-1 % searching the index of the accompanying field
                if rem(j,2)~=0 % fieldnames are placed on odd index
                    if strcmp(chklist{index},varargin{j})
                        I=j;
                    end
                end
            end
            options=setfield(options,chklist{index},varargin{I+1});
            index=[];
        end
        counter=counter+1;
    end
    options.k=options.k;
    options.sca=options.sca;
    options.loc=options.loc;  
end

if n==1 & p==1
    out=x;     %when X is a one by one matrix, all location estimators must be equal to that matrix
    return
elseif n==1
    x=x';      %we only want to work with column vectors
    n=p;
    p=1;
end

if n==2  % all location estimators must equal the average for n=2
    out=mean(x,1);
    return
end

alpha=0.413241928283814; %quadl('1/2.*sech(x/2).^2.*normpdf(x,0,1)',-10,10,1.e-15)
for i=1:p  
    X=x(:,i);
    t_0=feval(options.loc,X);
    s_0=feval(options.sca,X);
    tstep=t_0;
    if (s_0==0)
        out(i)=t_0;
    else
        j=1;
        while j<=options.k 
            z=(X-tstep)/s_0;
            y=tanh(z/2);
            tstep=tstep+s_0*(sum(y)/(n*alpha)); % updating location estimate
            y=[]; %clear help vector
            j=j+1;
        end
        out(i)=tstep;
    end
end

result = out;