l1median.m

prm,是一中对奇异样本不敏感的,偏最小二乘回归方法,已经有很多文章根据这个算法写出来了.使用时记得要引用chemlab上的文献.

function [mX]=l1median(X,tol);

% L1MEDIAN calculates the multivariate L1 median 
% I/O: [mX]=L1median(X,tol);
%
% X is the data matrix 
% tol is the convergence criterium; the iterative proces stops when ||m_k - m_{k+1}|| < tol.
%
% Ref: Hossjer and Croux (1995) "Generalizing Univariate Signed Rank Statistics for Testing
% and Estimating a Multivariate Location Parameter", Non-parametric Statistics, 4, 293-308.
%
% Written by C. Croux in GAUSS.
% Translated to MATLAB by S. Verboven.
% Adapted to Octave by S. Serneels, Uni Antwerp and TU Vienna, november 2003.

if nargin <2
   tol=1.e-02;
endif;
[n,p]=size(X);
maxstep=200;
%initializing starting value for m
m=median(X);
k=1;
while (k<=maxstep)
   mold=m;
   Xext=sortrows([norme(X-repmat(m,n,1)) X],1);
   dx=Xext(:,1);
   X=Xext(:,2:p+1);
   if all(dx)
      w=1./dx;
   else
      ww=dx(all(dx,2));
      w=1./ww;
      w=[zeros(length(dx)-length(w),1);w];
   endif
   delta=sum((X-repmat(m,n,1)).*repmat(w,1,p),1)./sum(w);
   nd=norme(delta);
   if all(nd<tol)
      maxhalf=0;
   else
      maxhalf=log2(nd/tol);
   endif
	m=mold+delta;   %computation of a new estimate
   nstep=0;
   while all(mrobj(X,m)>=mrobj(X,mold))&(nstep<=maxhalf)
      nstep=nstep+1;
      m=mold+delta./(2^nstep);
   endwhile
   if (nstep>maxhalf)
      mX=mold;
      break
   endif
   k=k+1;
endwhile
if k>maxstep
   disp('Iteration failed')
endif
mX=m;

endfunction



%-----
function n=norme(X);

%NORME calculates the euclidian norm of matrix X
% the output is a columnvector containing the norm of each row
%I/O: n=norme(X);

for i=1:size(X,1)
   n(i)=norm(X(i,:));
endfor;
n=n';

endfunction

%--------
function s=mrobj(X,m)

%MROBJ computes objective function in m based on X and a

xm=norme(X-repmat(m,size(X,1),1));
s=sum(xm,1)';

endfunction