irpls.m

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

%function [sam_weight,ir_residuals,Q,minfactor]=irpls(cal,val,caltar,valtar,iterations,factor)
function [sam_weight,ir_residuals,Q,minfactor]=irpls(cal,val,caltar,valtar,iterations,factor)
%To calculate the robust weights of each calibration sample

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
[ax,bx]=size(cal);
[vax,vbx]=size(val);
Qw=eye(ax);
oldcal=cal; oldcaltar=caltar;
oldval=val;
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
index=[];

for i=1:iterations
    Q=[];
    cal=Qw*oldcal(:,:);
    val=oldval;
    caltar=Qw*oldcaltar;
    [calrow,calcol]=size(cal);
    
    [res,valpress]=model_cen(cal,caltar);
    %factor=6;
    if nargin<6
       faco=valpress./min(valpress); 
       for ff=1:length(valpress)
            if faco(ff)<1.1
               factor=ff;
                break;
            end
        end
    end
    rmsecv(i)=sqrt(valpress(factor)/length(caltar));
    [pp,qq,ww,bb] = pls(cal', caltar', factor);
    [pred_tar, left_tar] = plspred(val', pp, qq, ww, bb,factor);
  
 rmsep=sqrt(sum((pred_tar-valtar').^2)/length(valtar));
 fprintf('rmsep=%f  \n',rmsep);
 %to calculate the weight of the objects
   result=res(:,factor);
   residuals=caltar-result;
   %media_resi=median(abs(residuals));
   media_resi=median(abs(residuals-median(residuals)));
    
   for ri=1:length(caltar)
       if abs(residuals(ri))<=6*media_resi
         Q(ri)=(1-(residuals(ri)/(6*media_resi)).^2).^2;  %Bisquate
        %  Q(ri)=1/(1+(residuals(ri)/(6*media_resi)).^2);   %Cauchy
        %Q(ri)=1/((1+abs(residuals(ri)/1.345*media_resi))).^2; %Fair
       %  Q(ri)=1;                                        %Huber     
      else  
          Q(ri)=0+10.^(-20);
       end
   end

   for ri=1:length(oldcaltar)
      if Q(ri)<0.001
         index=[index ri];
      end
   end
  Q(index)=0+10.^(-20);
  %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%5%%%%%%%%%%%%%%%%%%%%%%%%%%%%ȱһ
  Q=length(oldcaltar)*sqrt(Q)./sum(sqrt(Q));
  Qw=diag(Q);
end

%weight=Q./sum(Q);  
% sam_weight(1)=0;
% for i=1:length(Q)
%     sam_weight(i+1)=sum(weight(1:i));
% end
%  ir_residuals=residuals;
%  minfactor=factor;