npls.m

多维数据分析,有nPLS,PARAFAC,TURKER等

               L1=reshape(Factors(l_idx2(j,1):l_idx2(j,2)),dimx(j),Fac);
               Z=kr(L1,Z);
            end
            ZtZ=Z'*Z;
            ZtX=Z'*X';
            OldLoad=reshape(Factors(lidx(i,1):lidx(i,2)),DimX(i),Fac);
            L=pfls(ZtZ,ZtX,DimX(i),const(i),OldLoad,DoWeight,Weights);
            Factors(lidx(i,1):lidx(i,2))=L(:);
         end
         x=zeros(prod(DimX([1:ii-1 ii+1:ord])),DimX(ii));  % Rotate X so the current last mode is the first
         x(:)=X;
         X=x';
      end
   else
      for ii=ord:-1:1
         if ii==ord;
            i=1;
         else
            i=ii+1;
         end
         idd=[i+1:ord 1:i-1];
         l_idx2=lidx(idd,:);
         dimx=DimX(idd);
         if ~FixMode(i)
            L1=reshape(Factors(l_idx2(1,1):l_idx2(1,2)),dimx(1),Fac);
            if ord>2
               L2=reshape(Factors(l_idx2(2,1):l_idx2(2,2)),dimx(2),Fac);
               Z=kr(L2,L1);
            else
               Z = L1;
            end
            for j=3:ord-1
               L1=reshape(Factors(l_idx2(j,1):l_idx2(j,2)),dimx(j),Fac);
               Z=kr(L1,Z);
            end
            OldLoad=reshape(Factors(lidx(i,1):lidx(i,2)),DimX(i),Fac);
            L=pfls(Z,X,DimX(i),const(i),OldLoad,DoWeight,Weights);
            Factors(lidx(i,1):lidx(i,2))=L(:);
         end
         x=zeros(prod(DimX([1:ii-1 ii+1:ord])),DimX(ii));
         x(:)=X;
         X=x';
         x(:)=Weights;
         Weights=x';
      end
   end
   
   % EVALUATE SOFAR
   % Convert to new format
   clear ff,id1 = 0;
   for i = 1:length(DimX) 
      id2 = sum(DimX(1:i).*Fac);
      ff{i} = reshape(Factors(id1+1:id2),DimX(i),Fac);
      id1 = id2;
   end
   model=nmodel(ff);
   model = reshape(model,DimX(1),prod(DimX(2:end)));
   if MissMeth  % Missing values present
      connew=model(id);
      X(id)=model(id);
      errold=err;
      errX=X-model;
      if DoWeight==0
         err=sum(sum(errX(idmiss2).^2));
      else
         err=sum(sum((Weights(idmiss2).*errX(idmiss2)).^2));
      end
   else
      errold=err;
      if DoWeight==0
         err=sum(sum((X-model).^2));
      else
         err=sum(sum((Weights.*(X-model)).^2));
      end
   end
   
   if err<1000*eps, % Getting close to the machine uncertainty => stop
      disp(' WARNING')
      disp(' The misfit is approaching the machine uncertainty')
      disp(' If pure synthetic data is used this is OK, otherwise if the')
      disp(' data elements are very small it might be appropriate ')
      disp(' to multiply the whole array by a large number to increase')
      disp(' numerical stability. This will only change the solution ')
      disp(' by a scaling constant')
      f = 0;
   else
      f=abs((err-errold)/err);
      if f<crit % Convergence: then check that constraints are fulfilled
         if any(const==2)|any(const==3) % If nnls or unimodality imposed
            for i=1:ord % Extract the 
               if const(i)==2|const(i)==3 % If nnls or unimodality imposed
                  Loadd = Factors(sum(DimX(1:i-1))*Fac+1:sum(DimX(1:i))*Fac);
                  if any(Loadd<0)
                     ConstraintsNotRight=1;
                  else
                     ConstraintsNotRight=0;
                  end
               end
            end
         end
      end
   end
   
   if it/showfit-round(it/showfit)==0
      if showfit~=-1,
         ShowPhi=ShowPhi+1;
         if ShowPhi==ShowPhiWhen,
            ShowPhi=0;
            if showfit~=-1,
               disp(' '),
               disp('    Tuckers congruence coefficient'),
               % Convert to new format
               clear ff,id1 = 0;
               for i = 1:length(DimX) 
                  id2 = sum(DimX(1:i).*Fac);ff{i} = reshape(Factors(id1+1:id2),DimX(i),Fac);id1 = id2;
               end
               [phi,out]=ncosine(ff,ff);
               disp(phi),
               if MissMeth
                  fprintf(' Change in estim. missing values %12.10f',norm(connew-conold)/norm(conold));
                  disp(' ')
                  disp(' ')
               end
               disp(' Sum-of-Squares   Iterations  Explained')
               disp(' of residuals                 variation')
            end
         end
         if DoWeight==0
            PercentExpl=100*(1-err/SSX);
         else
            PercentExpl=100*(1-sum(sum((X-model).^2))/SSX);
         end
         fprintf(' %12.10f       %g        %3.4f    \n',err,it,PercentExpl);
         if Plt==2
            % Convert to new format
            clear ff,id1 = 0;
            for i = 1:length(DimX) 
               id2 = sum(DimX(1:i).*Fac);ff{i} = reshape(Factors(id1+1:id2),DimX(i),Fac);id1 = id2;
            end
            pfplot(reshape(X,DimX),ff,Weights',[0 0 0 0 0 0 0 1]);
            drawnow
         end
      end
   end
   
   
   
   % Make safety copy of loadings and initial parameters in temp.mat
   if it/50-round(it/50)==0
      save temp Factors
   end
   
   % JUDGE FIT
   if err>errold
      NumberOfInc=NumberOfInc+1;
   end
   
end % while f>crit


% CALCULATE TUCKERS CONGRUENCE COEFFICIENT
if showfit~=-1 & DimX(1)>1
   disp(' '),disp('   Tuckers congruence coefficient')
   % Convert to new format
   clear ff,id1 = 0;
   for i = 1:length(DimX) 
      id2 = sum(DimX(1:i).*Fac);ff{i} = reshape(Factors(id1+1:id2),DimX(i),Fac);id1 = id2;
   end
   [phi,out]=ncosine(ff,ff);
   disp(phi)
   disp(' ')
   if max(max(abs(phi)-diag(diag(phi))))>.85
      disp(' ')
      disp(' ')
      disp(' WARNING, SOME FACTORS ARE HIGHLY CORRELATED.')
      disp(' ')
      disp(' You could decrease the number of components. If this')
      disp(' does not help, try one of the following')
      disp(' ')
      disp(' - If systematic variation is still present you might')
      disp('   wanna decrease your convergence criterion and run')
      disp('   one more time using the loadings as initial guess.')
      disp(' ')
      disp(' - Or use another preprocessing (check for constant loadings)')
      disp(' ')
      disp(' - Otherwise try orthogonalising some modes,')
      disp(' ')
      disp(' - Or use Tucker3/Tucker2,')
      disp(' ')
      disp(' - Or a PARAFAC with some modes collapsed (if # modes > 3)')
      disp(' ')
   end
end


% SHOW FINAL OUTPUT

if DoWeight==0
   PercentExpl=100*(1-err/SSX);
else
   PercentExpl=100*(1-sum(sum((X-model).^2))/SSX);
end
if showfit~=-1
   fprintf(' %12.10f       %g        %3.4f \n',err,it,PercentExpl);
   if NumberOfInc>0
      disp([' There were ',num2str(NumberOfInc),' iterations that increased fit']);
   end
end


% POSTPROCES LOADINGS (ALL VARIANCE IN FIRST MODE)
A=reshape(Factors(lidx(1,1):lidx(1,2)),DimX(1),Fac);
for i=2:ord
   B=reshape(Factors(lidx(i,1):lidx(i,2)),DimX(i),Fac);
   for ff=1:Fac
      A(:,ff)=A(:,ff)*norm(B(:,ff));
      B(:,ff)=B(:,ff)/norm(B(:,ff));
   end
   Factors(lidx(i,1):lidx(i,2))=B(:);
end
Factors(lidx(1,1):lidx(1,2))=A(:);
if showfit~=-1
   disp(' ')
   disp(' Components have been normalized in all but the first mode')
end

% PERMUTE SO COMPONENTS ARE IN ORDER AFTER VARIANCE DESCRIBED (AS IN PCA) IF NO FIXED MODES
if ~any(FixMode)
   A=reshape(Factors(lidx(1,1):lidx(1,2)),DimX(1),Fac);
   [out,order]=sort(diag(A'*A));
   order=flipud(order);
   A=A(:,order);
   Factors(lidx(1,1):lidx(1,2))=A(:);
   for i=2:ord
      B=reshape(Factors(lidx(i,1):lidx(i,2)),DimX(i),Fac);
      B=B(:,order);
      Factors(lidx(i,1):lidx(i,2))=B(:);
   end  
   if showfit~=-1
      disp(' Components have been ordered according to contribution')
   end
elseif showfit ~= -1
   disp(' Some modes fixed hence no sorting of components performed')
end

% APPLY SIGN CONVENTION IF NO FIXED MODES


%  FixMode=1
if ~any(FixMode)&~(any(const==2)|any(const==3))
   Sign = ones(1,Fac);
   for i=ord:-1:2
      A=reshape(Factors(lidx(i,1):lidx(i,2)),DimX(i),Fac);
      Sign2=ones(1,Fac);
      for ff=1:Fac
         [out,sig]=max(abs(A(:,ff)));
         Sign(ff) = Sign(ff)*sign(A(sig,ff));
         Sign2(ff) = sign(A(sig,ff));
      end
      A=A*diag(Sign2);
      Factors(lidx(i,1):lidx(i,2))=A(:);
   end 
   A=reshape(Factors(lidx(1,1):lidx(1,2)),DimX(1),Fac);
   A=A*diag(Sign);
   Factors(lidx(1,1):lidx(1,2))=A(:);
   if showfit~=-1
      disp(' Components have been reflected according to convention')
   end
   
end 

% TOOLS FOR JUDGING SOLUTION
if nargout>3      
   x=X;
   if MissMeth
      x(id)=NaN*id;
   end
   % Convert to new format
   clear ff,id1 = 0;
   for i = 1:length(DimX) 
      id2 = sum(DimX(1:i).*Fac);ff{i} = reshape(Factors(id1+1:id2),DimX(i),Fac);id1 = id2;
   end
   corcondia=corcond(reshape(x,DimX),ff,Weights,1);
end

if Plt==1|Plt==2
   % Convert to new format
   clear ff,id1 = 0;
   for i = 1:length(DimX) 
      id2 = sum(DimX(1:i).*Fac);ff{i} = reshape(Factors(id1+1:id2),DimX(i),Fac);id1 = id2;
   end

   pfplot(reshape(X,DimX),ff,Weights,ones(1,8));
end

% Show which criterion stopped the algorithm
if showfit~=-1
   if ((f<crit) & (norm(connew-conold)/norm(conold)<MissConvCrit))
      disp(' The algorithm converged')
   elseif it==maxit
      disp(' The algorithm did not converge but stopped because the')
      disp(' maximum number of iterations was reached')
   elseif f<eps
      disp(' The algorithm stopped because the change in fit is now')
      disp(' smaller than the machine uncertainty.')
   else
      disp(' Algorithm stopped for some mysterious reason')
   end
end

% Convert to new format
clear ff,id1 = 0;
for i = 1:length(DimX) 
   id2 = sum(DimX(1:i).*Fac);ff{i} = reshape(Factors(id1+1:id2),DimX(i),Fac);id1 = id2;
end
Factors = ff;


function [A,B,C,fit]=dtld(X,F,SmallMode);

%DTLD direct trilinear decomposition
%
% See also:
% 'gram', 'parafac'
%
% Copyright, 1998 - 
% This M-file and the code in it belongs to the holder of the
% copyrights and is made public under the following constraints:
% It must not be changed or modified and code cannot be added.
% The file must be regarded as read-only. Furthermore, the
% code can not be made part of anything but the 'N-way Toolbox'.
% In case of doubt, contact the holder of the copyrights.
%
% Rasmus Bro
% Chemometrics Group, Food Technology
% Department of Food and Dairy Science
% Royal Veterinary and Agricultutal University
% Rolighedsvej 30, DK-1958 Frederiksberg, Denmark
% Phone  [phone]
% Fax    [phone]
% E-mail [e-mail]
%
%
% DIRECT TRILINEAR DECOMPOSITION
%
% calculate the parameters of the three-
% way PARAFAC model directly. The model
% is not the least-squares but will be close
% to for precise data with little model-error
%
% This implementation works with an optimal
% compression using least-squares Tucker3 fitting
% to generate two pseudo-observation matrices that
% maximally span the variation of all samples. per
% default the mode of smallest dimension is compressed
% to two samples, while the remaining modes are 
% compressed to dimension F.
% 
% For large arrays it is fastest to have the smallest
% dimension in the first mode
%
% INPUT
% [A,B,C]=dtld(X,F);
% X is the I x J x K array
% F is the number of factors to fit
% An optional parameter may be given to enforce which
% mode is to be compressed to dimension two
%
% Copyright 1998
% Rasmus Bro, KVL
% [e-mail]

% $ Version 2.00 $ May 2001 $ Changed to array notation $ RB $ Not compiled $
% $ Version 1.02 $ Date 28. July 1998 $ Not compiled $
% $ Version 1.03 $ Date 25. April 1999 $ Not compiled $

DimX = size(X);
X = reshape(X,DimX(1),prod(DimX(2:end)));

DontShowOutput = 1;

%rearrange X so smallest dimension is in first mode


if nargin<4
  [a,SmallMode] = min(DimX);
  X = nshape(reshape(X,DimX),SmallMode);
  DimX = DimX([SmallMode 1:SmallMode-1 SmallMode+1:3]);
  Fac   = [2 F F];
else
  X = nshape(reshape(X,DimX),SmallMode);
  DimX = DimX([SmallMode 1:SmallMode-1 SmallMode+1:3]);
  Fac   = [2 F F];
end
f=F;
if F==1;
  Fac   = [2 2 2];
  f=2;
end 


if DimX(1) < 2
  error(' The smallest dimension must be > 1')
end

if any(DimX(2:3)-Fac(2:3)<0)
  error(' This algorithm requires that two modes are of dimension not less the number of components')
end



% Compress data into a 2 x F x F array. Only 10 iterations are used since exact SL fit is insignificant; only obtaining good truncated bases is important
[Factors,Gt]=tucker(reshape(X,DimX),Fac,[0 0 0 0 NaN 10]);
% Convert to old format
Gt = reshape(Gt,size(Gt,1),prod(size(Gt))/size(Gt,1));

[At,Bt,Ct]=fac2let(Factors);

% Fit GRAM to compressed data
[Bg,Cg,Ag]=gram(reshape(Gt(1,:),f,f),reshape(Gt(2,:),f,f),F);

% De-compress data and find A


BB = Bt*Bg;
CC = Ct*Cg;
AA = X*pinv(kr(CC,BB)).';

if SmallMode == 1
  A=AA;
  B=BB;
  C=CC;
elseif SmallMode == 2 
  A=BB;
  B=AA;
  C=CC;
elseif SmallMode == 3
  A=BB;
  B=CC;
  C=AA;
end

fit = sum(sum(abs(X - AA*kr(CC,BB).').^2));
if ~DontShowOutput
  disp([' DTLD fitted raw data with a sum-squared error of ',num2str(fit)])
end


function mwa = outerm(facts,lo,vect)

if nargin < 2
  lo = 0;
end
if nargin < 3
  vect = 0;
end
order = length(facts);
if lo == 0
  mwasize = zeros(1,order);
else
  mwasize = zeros(1,order-1);
end
k = 0;
for i = 1:order
  if i ~= lo
    [m,n] = size(facts{i});
    k = k + 1;
    mwasize(k) = m;
    if k > 1
      if nofac ~= n
        error('All orders must have the same number of factors')
      end
    else
      nofac = n;
    end
  end
end
mwa = zeros(prod(mwasize),nofac);

for j = 1:nofac
  if lo ~= 1
    mwvect = facts{1}(:,j);
    for i = 2:order
	  if lo ~= i
        %mwvect = kron(facts{i}(:,j),mwvect);
		mwvect = mwvect*facts{i}(:,j)';
		mwvect = mwvect(:);
	  end
    end
  elseif lo == 1
    mwvect = facts{2}(:,j);
	for i = 3:order
      %mwvect = kron(facts{i}(:,j),mwvect);
	  mwvect = mwvect*facts{i}(:,j)';
	  mwvect = mwvect(:);
	end
  end
  mwa(:,j) = mwvect;
end
% If vect isn't one, sum up the results of the factors and reshape
if vect ~= 1
  mwa = sum(mwa,2);
  mwa = reshape(mwa,mwasize);
end