npls.m

多元统计分析是一种应用非常广泛的数据处理方法

function [Xfactors,Yfactors,Core,B,ypred,ssx,ssy,reg] = npls(X,Y,Fac,show);

%NPLS multilinear partial least squares regression
%
% See also:
% 'parafac' 'tucker'
%
%
% MULTILINEAR PLS  -  N-PLS
%
% INPUT
% X        Array of independent variables
% Y        Array of dependent variables
% Fac      Number of factors to compute
% 
% OPTIONAL
% show	   If show = NaN, no outputs are given
%
%
% OUTPUT
% Xfactors Holds the components of the model of X in a cell array.
%          Use fac2let to convert the parameters to scores and
%          weight matrices. I.e., for a three-way array do
%          [T,Wj,Wk]=fac2let(Xfactors);
% Yfactors Similar to Xfactors but for Y
% Core     Core array used for calculating the model of X
% B        The regression coefficients from which the scores in
%          the Y-space are estimated from the scores in the X-
%          space (U = TB);
% ypred    The predicted values of Y for one to Fac components
%          (array with dimension Fac in the last mode)
% ssx      Variation explained in the X-space.
%          ssx(f+1,1) is the sum-squared residual after first f factors.
%          ssx(f+1,2) is the percentage explained by first f factors.
% ssy      As above for the Y-space
% reg      Cell array with regression coefficients for raw (preprocessed) X
%
%
% AUXILIARY
%
% If missing elements occur these must be represented by NaN.
%
%
% [Xfactors,Yfactors,Core,B,ypred,ssx,ssy,reg] = npls(X,y,Fac);
% or short
% [Xfactors,Yfactors,Core,B] = npls(X,y,Fac);
%

% Copyright (C) 1995-2006  Rasmus Bro & Claus Andersson
% Copenhagen University, DK-1958 Frederiksberg, Denmark, [e-mail]
%
% This program is free software; you can redistribute it and/or modify it under 
% the terms of the GNU General Public License as published by the Free Software 
% Foundation; either version 2 of the License, or (at your option) any later version.
%
% This program is distributed in the hope that it will be useful, but WITHOUT 
% ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or FITNESS 
% FOR A PARTICULAR PURPOSE.  See the GNU General Public License for more details.
% You should have received a copy of the GNU General Public License along with 
% this program; if not, write to the Free Software Foundation, Inc., 51 Franklin 
% Street, Fifth Floor, Boston, MA  02110-1301, USA.


% $ Version 1.02 $ Date July 1998 $ Not compiled $
% $ Version 1.03 $ Date 4. December 1998 $ Not compiled $ Cosmetic changes
% $ Version 1.04 $ Date 4. December 1999 $ Not compiled $ Cosmetic changes
% $ Version 1.05 $ Date July 2000 $ Not compiled $ error caused weights not to be normalized for four-way and higher
% $ Version 1.06 $ Date November 2000 $ Not compiled $ increase max it and decrease conv crit to better handle difficult data
% $ Version 2.00 $ May 2001 $ Changed to array notation $ RB $ Not compiled $
% $ Version 2.01 $ June 2001 $ Changed to handle new core in X $ RB $ Not compiled $
% $ Version 2.02 $ January 2002 $ Outputs all predictions (1 - LV components) $ RB $ Not compiled $
% $ Version 2.03 $ March 2004 $ Changed initialization of u $ RB $ Not compiled $
% $ Version 2.04 $ Jan 2005 $ Modified sign conventions of scores and loads $ RB $ Not compiled $
% $ Version 3.00 $ Aug 2007 $ Serious error in sign switch fixed $ RB $ Not compiled $

if nargin==0
   disp(' ')
   disp(' ')
   disp(' THE N-PLS REGRESSION MODEL')
   disp(' ')
   disp(' Type <<help npls>> for more info')
   disp('  ')
   disp(' [Xfactors,Yfactors,Core,B,ypred,ssx,ssy] = npls(X,y,Fac);')
   disp(' or short')
   disp(' [Xfactors,Yfactors,Core,B] = npls(X,y,Fac);')
   disp(' ')
   return
elseif nargin<3
   error(' The inputs X, y, and Fac must be given')
end



if ~exist('show')==1|nargin<4
   show=1;
end

maxit=120;

DimX = size(X);
X = reshape(X,DimX(1),prod(DimX(2:end)));
ordX = length(DimX);if ordX==2&size(X,2)==1;ordX = 1;end
DimY = size(Y);
Y = reshape(Y,DimY(1),prod(DimY(2:end)));
ordY = length(DimY);if ordY==2&size(Y,2)==1;ordY = 1;end


[I,Jx]=size(X);
[I,Jy]=size(Y);

missX=0;
missy=0;
MissingX = 0;
MissingY = 0;
if any(isnan(X(:)))|any(isnan(Y(:)))
   if any(isnan(X(:)))
      MissingX=1;
   else
      MissingX=0;
   end
   if any(isnan(Y(:)))
      MissingY=1;
   else
      MissingY=0;
   end
   if show~=0&~isnan(show)
      disp(' ')
      disp(' Don''t worry, missing values will be taken care of')
      disp(' ')
   end
   missX=abs(1-isnan(X));
   missy=abs(1-isnan(Y));
end
crit=1e-10;
B=zeros(Fac,Fac);
T=[];
U=[];
Qkron =[];
if MissingX
   SSX=sum(sum(X(find(missX)).^2));
else
   SSX=sum(sum(X.^2));
end
if MissingY
   SSy=sum(sum(Y(find(missy)).^2));
else
   SSy=sum(sum(Y.^2));
end
ssx=[];
ssy=[];
Xres=X;
Yres=Y;
xmodel=zeros(size(X));
Q=[];
W=[];

for num_lv=1:Fac
   
   %init
   
   
   % u=rand(DimX(1),1); Old version
   if size(Yres,2)==1
     u = Yres;
   else
     [u] = pcanipals(Yres,1,0);   
   end
   
   t=rand(DimX(1),1);
   tgl=t+2;it=0;
   while (norm(t-tgl)/norm(t))>crit&it<maxit
      tgl=t;
      it=it+1;
      
      % w=X'u
      [wloads,wkron] = Xtu(X,u,MissingX,missX,Jx,DimX,ordX);
      
      % t=Xw
      if MissingX
         for i=1:I,
            m = find(missX(i,:));
            t(i)=X(i,m)*wkron(m)/(wkron(m)'*wkron(m));
         end
      else
         t=X*wkron;
      end
      
      % w=X'u
      [qloads,qkron] = Xtu(Yres,t,MissingY,missy,Jy,DimY,ordY);
      % u=yq
      if MissingY
         for i=1:I
            m = find(missy(i,:));
            u(i)=Yres(i,m)*qkron(m)/(qkron(m)'*qkron(m));
         end
      else
         u=Yres*qkron;
      end
   end
   
   
%    % Fix signs
%    [Factors] = signswtch({t,wloads{:}},X);
%    t = Factors{1};
%    wloads = Factors(2:end);
%    % Fix signs
%    [Factors] = signswtch({u,qloads{:}},X);
%    u = Factors{1};
%    qloads = Factors(2:end);
   
   % Arrange t scores so they positively correlated with u
   cc = corrcoef([t u]);
   if sign(cc(2,1))<0
     t = -t;
     wloads{1}=-wloads{1};
   end

   
   
   
   T=[T t];
   for i = 1:ordX-1
      if num_lv == 1
         W{i} = wloads{i};
      else
         W{i} = [W{i} wloads{i}];
      end
   end
   U=[U u];
   for i = 1:max(ordY-1,1)
      if num_lv == 1
         Q{i} = qloads{i};
      else
         Q{i} = [Q{i} qloads{i}];
      end
   end
   Qkron = [Qkron qkron];
  
   % Make core arrays
   if ordX>1
      Xfac{1}=T;Xfac(2:ordX)=W;
      Core{num_lv} = calcore(reshape(X,DimX),Xfac,[],0,1);
   else
      Core{num_lv} = 1;
   end
%   if ordY>1
%      Yfac{1}=U;Yfac(2:ordY)=Q;
%      Ycore{num_lv} = calcore(reshape(Y,DimY),Yfac,[],0,1);
%   else
%      Ycore{num_lv} = 1;
%   end
   
   
   B(1:num_lv,num_lv)=inv(T'*T)*T'*U(:,num_lv);
   
   if Jy > 1
      if show~=0&~isnan(show)
         disp(' ') 
         fprintf('number of iterations: %g',it);
         disp(' ')
      end
   end
   
   % Make X model
   if ordX>2
      Wkron = kron(W{end},W{end-1});
   else
      Wkron = W{end};
   end
   for i = ordX-3:-1:1
      Wkron = kron(Wkron,W{i});
   end
   if num_lv>1
      xmodel=T*reshape(Core{num_lv},num_lv,num_lv^(ordX-1))*Wkron';
   else
      xmodel = T*Core{num_lv}*Wkron';
   end
   
   % Make Y model   
 %  if ordY>2
 %     Qkron = kron(Q{end},Q{end-1});
 %  else
 %     Qkron = Q{end};
 %  end
 %  for i = ordY-3:-1:1
 %     Qkron = kron(Qkron,Q{i});
 %  end
 %  if num_lv>1
 %     ypred=T*B(1:num_lv,1:num_lv)*reshape(Ycore{num_lv},num_lv,num_lv^(ordY-1))*Qkron';
 %  else
 %     ypred = T*B(1:num_lv,1:num_lv)*Ycore{num_lv}*Qkron';
 %  end
 ypred=T*B(1:num_lv,1:num_lv)*Qkron';
 Ypred(:,num_lv) = ypred(:); % Vectorize to avoid problems with different orders and the de-vectorize later on
   
   Xres=X-xmodel; 
   Yres=Y-ypred;
   if MissingX
      ssx=[ssx;sum(sum(Xres(find(missX)).^2))];
   else
      ssx=[ssx;sum(sum(Xres.^2))];
   end
   if MissingY
      ssy=[ssy;sum(sum((Y(find(missy))-ypred(find(missy))).^2))];
   else
      ssy=[ssy;sum(sum((Y-ypred).^2))];
   end
end
ypred = reshape(Ypred',[size(Ypred,2) DimY]);
ypred = permute(ypred,[2:ordY+1 1]);

ssx= [ [SSX(1);ssx] [0;100*(1-ssx/SSX(1))]];
ssy= [ [SSy(1);ssy] [0;100*(1-ssy/SSy(1))]];

if show~=0&~isnan(show)
   disp('  ')
   disp('   Percent Variation Captured by N-PLS Model   ')
   disp('  ')
   disp('   LV      X-Block    Y-Block')
   disp('   ----    -------    -------')
   ssq = [(1:Fac)' ssx(2:Fac+1,2) ssy(2:Fac+1,2)];
   format = '   %3.0f     %6.2f     %6.2f';
   for i = 1:Fac
      tab = sprintf(format,ssq(i,:)); disp(tab)
   end
end

Xfactors{1}=T;
for j = 1:ordX-1
   Xfactors{j+1}=W{j};
end

Yfactors{1}=U;
for j = 1:max(ordY-1,1)
   Yfactors{j+1}=Q{j};
end


% Calculate regression coefficients that apply directly to X
  if nargout>7
    if length(DimY)>2
      error(' Regression coefficients are only calculated for models with vector Y or multivariate Y (not multi-way Y)')
    end
      R = outerm(W,0,1);
    for iy=1:size(Y,2)
      if length(DimX) == 2
        dd = [DimX(2) 1];
      else
        dd = DimX(2:end);
      end
      for i=1:Fac
        sR = R(:,1:i)*B(1:i,1:i)*diag(Q{1}(iy,1:i));
        ssR = sum( sR',1)';
        reg{iy,i} = reshape( ssR ,dd);
      end       
    end
    
  end




function [wloads,wkron] = Xtu(X,u,Missing,miss,J,DimX,ord);


% w=X'u
if Missing
   for i=1:J
      m = find(miss(:,i));
      if (u(m)'*u(m))~=0
        ww=X(m,i)'*u(m)/(u(m)'*u(m));
      else
        ww=X(m,i)'*u(m);
      end
      if length(ww)==0
         w(i)=0;
      else
         w(i)=ww;
      end
   end
else
   w=X'*u;
end

% Reshape to array
if length(DimX)>2
   w_reshaped=reshape(w,DimX(2),prod(DimX(3:length(DimX))));
else
   w_reshaped = w(:);
end


% Find one-comp decomposition
if length(DimX)==2
   wloads{1} = w_reshaped/norm(w_reshaped);
   
elseif length(DimX)==3&~any(isnan(w_reshaped))
   [w1,s,w2]=svd(w_reshaped);
   wloads{1}=w1(:,1);
   wloads{2}=w2(:,1);
else
   wloads=parafac(reshape(w_reshaped,DimX(2:length(DimX))),1,[0 2 0 0 NaN]');
   for j = 1:length(wloads);
      wloads{j} = wloads{j}/norm(wloads{j});
   end
end

% Apply sign convention
for i = 1:length(wloads)
   sq = (wloads{i}.^2).*sign(wloads{i});
   wloads{i} = wloads{i}*sign(sum(sq));