一些偏最小二乘法(PLS)间隔偏最小二乘(siPLS)联合区间偏最小二乘等的建模代码。

源代码在线查看: pls_val.m

软件大小: 1146 K
上传用户: zhangtaoai007
关键词: siPLS PLS 最小二乘法 建模
下载地址: 免注册下载 普通下载 VIP

相关代码

				function PLSmodel = pls_val(X,Y,no_of_lv,prepro_method,val_method,segments);
				
				%  pls_val for PLS modelling with selected validation method
				%
				%  Input:
				%  X contains the independent variables
				%  Y contains the dependent variable(s), NOTE: Y is allways autoscaled
				%  no_of_lv is the number of PLS components
				%  prepro_method is 'mean', 'auto', 'mscmean', 'mscauto' or 'none'
				%  val_method is 'test', 'full', 'syst111', 'syst123', 'random' or 'manual'
				%  segments is number of segments in cross validation
				%       if val_method is 'test' then segments should be a column vector with test set indices
				%       if val_method is 'manual' then segments should be a cell array, see makeManualSegments
				%  Output:
				%  PLSmodel is a structured array containing all model and validation information
				%
				%  Subfunctions at the end of this file: rmbi, msc, msc_pre, mncn, auto
				%
				%  Copyright, Chemometrics Group - KVL, Copenhagen, Denmark
				%  Lars N鴕gaard, July 2004
				%
				%  PLSmodel = pls_val(X,Y,no_of_lv,prepro_method,val_method,segments);
				
				[n,m] = size(X);
				[o,p] = size(Y);
				
				if nargin==6 & max(size(segments))==1 % Not test or manual
				  no_sampl=fix(n/segments);
				  left_over_samples=mod(n,segments);
				end
				
				if nargin==6 & iscell(segments)
				  manualseg=segments; % A cell array
				  segments=max(size(segments)); % Now a scalar
				  Nsamples=0;
				  for j=1:segments
				    Nsamples=Nsamples + max(size(manualseg{j}));
				  end
				  if n~=Nsamples
				    disp('The number of samples in X does not correspond to the number of samples in manualseg')
				    return
				  end
				end
				
				% if strcmp(lower(val_method),'full')
				%     if nargin==6 & ~isempty(segments)
				%         disp('No need to include segments when full cross validation is chosen')
				%     end
				% end
				
				Ypred=zeros(n,p,no_of_lv);
				count=1;
				
				PLSmodel.prepro_method=prepro_method;
				PLSmodel.val_method=val_method;
				
				if strcmp(val_method,'full')
				   val_method='syst123';
				   segments=n;
				end
				
				if strcmp(val_method,'random')
				    ix=rand(n,1);
				    [temp,ix]=sort(ix);
				end
				
				switch val_method
				    case 'test'
				      tot=(1:n)';
				      tot(segments)=[];
				      cal = tot;
				      Xseg=X(cal,:);
				      Yseg=Y(cal,:);
				      Xpseg=X(segments,:);
				      [Yseg,my,stdy]=auto(Yseg);
				      if strcmp(lower(prepro_method),'mean')
				           [Xseg,mx]=mncn(Xseg);
				           Xpseg=scalenew(Xpseg,mx);
				      elseif strcmp(lower(prepro_method),'auto')
				           [Xseg,mx,stdx]=auto(Xseg);
				           Xpseg=scalenew(Xpseg,mx,stdx);
				      elseif strcmp(lower(prepro_method),'mscmean')
				           [Xseg,Xsegmeancal]=msc(Xseg);
				           [Xseg,mx]=mncn(Xseg);
				           Xpseg=msc_pre(Xpseg,Xsegmeancal);
				           Xpseg=scalenew(Xpseg,mx);
				      elseif strcmp(lower(prepro_method),'mscauto')
				           [Xseg,Xsegmeancal]=msc(Xseg);
				           [Xseg,mx,stdx]=auto(Xseg);
				           Xpseg=msc_pre(Xpseg,Xsegmeancal);
				           Xpseg=scalenew(Xpseg,mx,stdx);
				      elseif strcmp(lower(prepro_method),'none')
				           % To secure that no centering/scaling is OK    
				      end
				      [P,Q,W,T,U,bsco,ssqdif]=pls(Xseg,Yseg,no_of_lv);
				      Ypred(segments,:,:)=pls_pre(Xpseg,bsco,P,Q,W,no_of_lv);
				      for j=1:no_of_lv
				        Ypred(segments,:,j)=scaleback(Ypred(segments,:,j),my,stdy);
				      end
				      Ypred0(segments,:)=ones(max(size(segments)),1)*zeros(size(my)); % For zero PLSC estimate as average of calibration segment
				      Ypred0(segments,:)=scaleback(Ypred0(segments,:),my,stdy);
				      
				      PLSmodel.test=segments;
				      
				    case {'full','syst111','syst123','random','manual'}
				    % SegmentBar = waitbar(0,'Working on cross-validation...','Name','Cross-validation...'); % NEW
					for i=1:segments
				        % waitbar(i/segments,SegmentBar,['Performing cross-validation on segment',num2str(i),' of ',num2str(segments)]); % NEW
				        if strcmp(val_method,'syst111')
				              if left_over_samples==0
				                   count=count;
				                   p_cvs=((i-1)*no_sampl+1+(count-1):i*no_sampl+(count-1))';
				              else   
				                   p_cvs=((i-1)*no_sampl+1+(count-1):i*no_sampl+count)';
				                   count=count+1;
				                   left_over_samples=left_over_samples-1;
				              end
				         elseif strcmp(val_method,'syst123')
				              p_cvs=(i:segments:n)';
				         elseif strcmp(val_method,'random')
				              p_cvs=(i:segments:n)';
				              p_cvs=ix(p_cvs)';
				         elseif strcmp(val_method,'manual')
				              p_cvs=manualseg{i};
				         end
				         tot=(1:n)';
				         tot(p_cvs)=[];
				         m_cvs = tot;
				         PLSmodel.cv{i}=p_cvs;
				         Xseg=X(m_cvs,:);
				         Yseg=Y(m_cvs,:);
				         Xpseg=X(p_cvs,:);
				         [Yseg,my,stdy]=auto(Yseg);
				         if strcmp(lower(prepro_method),'mean')
				              [Xseg,mx]=mncn(Xseg);
				              Xpseg=scalenew(Xpseg,mx);
				         elseif strcmp(lower(prepro_method),'auto')
				              [Xseg,mx,stdx]=auto(Xseg);
				              Xpseg=scalenew(Xpseg,mx,stdx);
				         elseif strcmp(lower(prepro_method),'mscmean')
				              [Xseg,Xsegmeancal]=msc(Xseg);
				              [Xseg,mx]=mncn(Xseg);
				              Xpseg=msc_pre(Xpseg,Xsegmeancal);
				              Xpseg=scalenew(Xpseg,mx);
				         elseif strcmp(lower(prepro_method),'mscauto')
				              [Xseg,Xsegmeancal]=msc(Xseg);
				              [Xseg,mx,stdx]=auto(Xseg);
				              Xpseg=msc_pre(Xpseg,Xsegmeancal);
				              Xpseg=scalenew(Xpseg,mx,stdx);
				         elseif strcmp(lower(prepro_method),'none')
				              % To secure that no centering/scaling is OK    
				         end
				         [P,Q,W,T,U,bsco,ssqdif]=pls(Xseg,Yseg,no_of_lv);
				         Ypred(p_cvs,:,:)=pls_pre(Xpseg,bsco,P,Q,W,no_of_lv);
				         for j=1:no_of_lv
				           Ypred(p_cvs,:,j)=scaleback(Ypred(p_cvs,:,j),my,stdy);
				         end
				         Ypred0(p_cvs,:)=ones(max(size(p_cvs)),1)*zeros(size(my)); % For zero PLSC estimate as average of calibration segment
				         Ypred0(p_cvs,:)=scaleback(Ypred0(p_cvs,:),my,stdy);
					end
					% delete(SegmentBar); % NEW
				end % switch val_method
				
				switch val_method
				  case 'test'
				    [RMSE0,Bias0]=rmbi(Y(segments,:),Ypred0(segments,:));
				  otherwise
				    [RMSE0,Bias0]=rmbi(Y,Ypred0);
				end
				
				for i=1:no_of_lv
				  switch val_method
				    case 'test'
				      [RMSE(i),Bias(i)]=rmbi(Y(segments,:),Ypred(segments,:,i)); % Subfunction at the end of this file
				    otherwise
				      [RMSE(i),Bias(i)]=rmbi(Y,Ypred(:,:,i)); % Subfunction at the end of this file
				  end
				end
				
				PLSmodel.Ypred0=Ypred0;
				PLSmodel.Ypred=Ypred;
				PLSmodel.RMSE=[RMSE0 RMSE];
				PLSmodel.Bias=[Bias0 Bias];
				
				% Global model with all samples
				[Y,my,stdy]=auto(Y);
				if strcmp(lower(prepro_method),'mean')
					[X,mx]=mncn(X);
				elseif strcmp(lower(prepro_method),'auto')
					[X,mx,stdx]=auto(X);
				elseif strcmp(lower(prepro_method),'mscmean')
					X=msc(X);
					[X,mx]=mncn(X);
				elseif strcmp(lower(prepro_method),'mscauto')
					X=msc(X);
					[X,mx,stdx]=auto(X);
				elseif strcmp(lower(prepro_method),'none')
				    disp('No scaling')
				    % To secure that no centering/scaling is OK    
				end
				
				[PLSmodel.P,PLSmodel.Q,PLSmodel.W,PLSmodel.T,PLSmodel.U,PLSmodel.bsco,PLSmodel.ssqdif]=pls(X,Y,no_of_lv);
				
				function [RMSE,Bias]=rmbi(Yref,Ypred)
				[n,m]=size(Yref);
				RMSE = sqrt( sum(sum((Ypred-Yref).^2))/(n*m) );
				Bias = sum(sum(Ypred-Yref))/(n*m);
				
				function [Xmsc,Xmeancal]=msc(X)
				[n,m]=size(X);
				Xmeancal=mean(X);
				for i=1:n
				  coef=polyfit(Xmeancal,X(i,:),1);
				  Xmsc(i,:)=(X(i,:)-coef(2))/coef(1);
				end
				
				function Xpmsc=msc_pre(Xp,Xmeancal)
				[n,m]=size(Xp);
				for i=1:n
				  coef=polyfit(Xmeancal,Xp(i,:),1);
				  Xpmsc(i,:)=(Xp(i,:)-coef(2))/coef(1);
				end
				
				function [Xmean,meanX] = mncn(X)
				[n,m] = size(X);
				meanX = mean(X);
				Xmean = (X-meanX(ones(n,1),:));
				
				function [Xauto,meanX,stdX] = auto(X)
				[n,m] = size(X);
				meanX = mean(X);
				stdX  = std(X);
				Xauto = (X-meanX(ones(n,1),:))./stdX(ones(n,1),:);
				
				function Xscalenew = scalenew(Xnew,meanXold,stdXold)
				[n,m] = size(Xnew);
				if nargin == 2
				  Xscalenew = (Xnew-meanXold(ones(n,1),:));
				elseif nargin == 3
				  Xscalenew = (Xnew-meanXold(ones(n,1),:))./stdXold(ones(n,1),:);
				end
				
				function Xscaleback = scaleback(X,meanX,stdX)
				[n,m] = size(X);
				if nargin == 2
				  Xscaleback = X + meanX(ones(n,1),:);
				elseif nargin == 3
				  Xscaleback = (X.*stdX(ones(n,1),:)) + meanX(ones(n,1),:);
				end
							

相关资源