prm.m

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

function [b,wy,wt,T,P,R]=prm(X,y,h,fairct,opt)

% PRM Partial Robust M-regression estimator
% ------------------------------------------
% Inputs:  y, vector of size (n,1): contains responses
%          X, matrix of size (n,p): contains predictor variables in columns
%          h, scalar:  the number of components
%          fairct,  scalar: the tuning constant of the weighting function (default=4)
%          opt, string:  centering done by L1-median if opt2="l1m", or  by
%                        coordinatewise median if opt2="cm" (default="l1m")
% ------------------------------------------
% Outputs: b, vector of size (p,1): contains  regression coefficients
%          wy, vector of size (n,1): contains residual weights
%          wt, vector of size (n,1): contains weights for leverage points
%          T, vector of size (n,h), contains the scores of the observations in the rows
%          P, matrix of size (p,h), contains the loading vectors in its colums
%          R, matrix of size (p,h), contains the PRM weighting vectors as 
%                                   columns
% ------------------------------------------
% The partial robust M-regression estimator is proposed in S. Serneels, C. Croux, P. Filzmoser and P.J. Van Espen,
% "Partial M-regression", submitted for publication.

% Written by Sven Serneels, University of Antwerp, and C. Croux, KULeuven.

%

% Uses the program Unisimpls.m
if nargin < 5
    opt='l1m';
    if nargin < 4
        fairct=4;
    end;
end;


[n,p]=size(X);

if p>n
    dimensions=1;
    dimension=p-n;
    [V,S,U]=svd(X',0);
    X=U*S;
    [n,p]=size(X);
else
    dimensions=0;
end;

if (opt=='l1m')
    mx=l1median(X,1e-5);
else
    mx=median(X);
end
my=median(y);

Xmc=X-repmat(mx,n,1);
wx=sqrt(sum(Xmc'.^2))';
wx=wx./median(wx);
wx=1./(1+abs(wx./fairct)).^2;

wy=abs(y-repmat(my,n,1));
wy=wy./median(wy);
wy=1./(1+abs(wy./fairct)).^2;

w=wx.*wy;
Xw=X.*repmat(sqrt(w),1,p);
yw=y.*sqrt(w);

loops=1;
ngamma=10^5;
difference=1;

while ~(any([difference < 10^-2 ; loops>30]))
    ngammaold=ngamma;
    [b,T,P,R]=Unisimpls(Xw,yw,Xw'*yw,h);
    b=b(:,h);
    gamma=(yw'*T)';
    T=T.*repmat(1./sqrt(w),1,h);
    r=y-T*gamma;
    rc=r-median(r);
    r=rc./median(abs(rc));
    wy=1./(1+abs(r./fairct)).^2;
    if opt=='l1m'
        mt=l1median(T,1e-5);
    else
        mt=median(T);
    end
    dt=T-repmat(mt,n,1);
    wt=sqrt(sum(dt.^2,2));
    wt=wt./median(wt);
    wt=1./(1+abs(wt./fairct)).^2;
    ngamma=sqrt(sum(gamma.^2));
    difference=abs(ngamma-ngammaold)./ngamma;
    w=wy.*wt;
    Xw=X.*repmat(sqrt(w),1,p);
    yw=y.*sqrt(w);
    loops=loops+1;
end;


if dimensions
    b=V*b;
end;
%
endfunction
% -------------------------------------------------------------------------------------------------
function [B,T,P,R,V] = unisimpls(X,y,s,A)

% function [T,P,R] = UNISIMPLS(X,y,s,A)
% Full implementation of SIMPLS algorithm for univariate y
% X and y are assumed to be centered
% s = X' * y
% A = # dimensions
% T is the X-block score matrix; the columns of T have unit length
% P = X' * T is the X-block loading matrix
% R is the matrix of simPLS weights

% Sijmen de Jong
% Unilever Research Laboratorium Vlaardingen

[n,px] = size(X);
T = zeros(n,A); R = zeros(px,A); P = R; V = R; v = zeros(px,1); B=V;
for a = 1:A
   r = s;
   t = X * r; t =  t - T(:,1:max(1,a-1)) * (T(:,1:max(1,a-1))' * t);
   normt =  sqrt(t'*t); t = t / normt; r = r / normt;
   p = X' * t;
   v = p - V(:,1:max(1,a-1)) * (V(:,1:max(1,a-1))' * p);
   v = v / sqrt(v'*v);
   s = s - v * (v' * s);
   T(:,a) = t; P(:,a) = p; R(:,a) = r; V(:,a) = v;
   B(:,a)=R(:,1:a)*R(:,1:a)'*X'*y;
end
clear t p r v;
endfunction