impplot.m

Jackknife PARAFAC, 用在多维线性分解.

    if DoWeight==0
      err=sum(sum((X-model).^2));
    else
      err=sum(sum((Weights.*(X-model)).^2));
    end
  end
  if err/SSX<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|Plt==3
        % 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
     % POSTPROCESS. IF PCA on two-way enforce orth in both modes.
   
end % while f>crit

   if DoingPCA 
       A=reshape(Factors(lidx(1,1):lidx(1,2)),DimX(1),Fac);
       B=reshape(Factors(lidx(i,1):lidx(i,2)),DimX(i),Fac);
       [u,s,v]=svd(A*B',0);
       A = u(:,1:size(A,2))*s(1:size(A,2),1:size(A,2));
       B = u(:,1:size(B,2));
       Factors = [A(:);B(:)];
   end


% 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)
if Options(4)==0|Options(4)==1
  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
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,0);
end

if Plt==1|Plt==2|Plt==3
  % 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
  %   if Fac<6&Plt~=3&order>2&ord>2
  if Fac<6&Plt~=3&ord>2
    pfplot(reshape(X,DimX),ff,Weights,ones(1,8));
  else
    pfplot(reshape(X,DimX),ff,Weights,[1 1 0 1 1 1 1 1]);
    if ord>2
      disp(' Core consistency plot not shown because it requires large memory')
      disp(' It can be made writing pfplot(X,Factors,[Weights],[0 0 1 0 0 0 0 0]');
    else
      disp(' Core consistency not applicable for two-way data')
    end
  end
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 AB = kr(A,B);
%KR Khatri-Rao product
%
% The Khatri - Rao product
% For two matrices with similar column dimension the khatri-Rao product
% is kr(A,B) = [kron(B(:,1),A(:,1) .... kron(B(:,F),A(:,F)]
% 
% I/O AB = ppp(A,B);
%
% kr(A,B) equals ppp(B,A) - where ppp is the triple-P product = 
% the parallel proportional profiles product which was originally 
% suggested in Bro, Ph.D. thesism, 1998

% 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]
%
% $ Version 1.02 $ Date 28. July 1998 $ Not compiled $
% $ Version 2.00 $ May 2001 $ Changed to array notation $ RB $ Not compiled $

[I,F]=size(A);
[J,F1]=size(B);

if F~=F1
   error(' Error in kr.m - The matrices must have the same number of columns')
end

AB=zeros(I*J,F);
for f=1:F
   ab=B(:,f)*A(:,f).';
   AB(:,f)=ab(:);
end


function load=pfls(ZtZ,ZtX,dimX,cons,OldLoad,DoWeight,W);

%PFLS
%
% See also:
% 'unimodal' 'monreg' 'fastnnls'
%
% 
% Calculate the least squares estimate of
% load in the model X=load*Z' => X' = Z*load'
% given ZtZ and ZtX
% cons defines if an unconstrained solution is estimated (0)
% or an orthogonal (1), a nonnegativity (2), or a unimodality (3)
%
%
% Used by PARAFAC.M

% $ Version 1.02 $ Date 28. July 1998 $ Not compiled $
%
% 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]
%

% Apr 2002 - Fixed error in weighted ls $ rb

if ~DoWeight

  if cons==0 % No constr
    %load=((Z'*Z)\Z'*Xinuse)';
    load=(pinv(ZtZ)*ZtX)';
  
  elseif cons==1 % Orthogonal loadings acc. to Harshman & Lundy 94
    load=ZtX'*(ZtX*ZtX')^(-.5);

  elseif cons==2 % Nonnegativity constraint
    load=zeros(size(OldLoad));
    for i=1:dimX
       load(i,:)=fastnnls(ZtZ,ZtX(:,i))';
%       if min(load(i,:))<-eps*1000
%          load(i,:)=OldLoad(i,:);
%       end
    end

  elseif cons==3 % Unimodality & NNLS
     load=OldLoad;
     F=size(OldLoad,2);
     if F>1
       for i=1:F
        ztz=ZtZ(i,i);
        ztX=ZtX(i,:)-ZtZ(i,[1:i-1 i+1:F])*load(:,[1:i-1 i+1:F])';
        beta=(pinv(ztz)*ztX)';
        load(:,i)=ulsr(beta,1);
       end
     else
       beta=(pinv(ZtZ)*ZtX)';
       load=ulsr(beta,1);
     end
  end

elseif DoWeight
  Z=ZtZ;
  X=ZtX;
  if cons==0 % No constr
    load=OldLoad;
    one=ones(1,size(Z,2));
    for i=1:dimX
      ZW=Z.*(W(i,:).^2'*one);
      %load(i,:)=(pinv(Z'*diag(W(i,:))*Z)*(Z'*diag(W(i,:))*X(i,:)'))';
      load(i,:)=(pinv(ZW'*Z)*(ZW'*X(i,:)'))';
    end

  elseif cons==2 % Nonnegativity constraint
    load=OldLoad;
    one=ones(1,size(Z,2));
    for i=1:dimX
      ZW=Z.*(W(i,:).^2'*one);
      load(i,:)=fastnnls(ZW'*Z,ZW'*X(i,:)')';
    end

  elseif cons==1
    disp(' Weighted orthogonality not implemented yet')
    disp(' Please contact the authors for further information')
    error

  elseif cons==3
    disp(' Weighted unimodality not implemented yet')
    disp(' Please contact the authors for further information')
    error

  end

end


% Check that NNLS and alike do not intermediately produce columns of only zeros
if cons==2|cons==3
  if any(sum(load)==0)  % If a column becomes only zeros the algorithm gets instable, hence the estimate is weighted with the prior estimate. This should circumvent numerical problems during the iterations
    load = .9*load+.1*OldLoad;
  end
end

function [Xm]=nmodel(Factors,G,Om);

%NMODEL make model of data from loadings
%
% function [Xm]=nmodel(Factors,G,Om);
%
% This algorithm requires access to:
% 'neye.m'
%
%
% [Xm]=nmodel(Factors,G,Om);
%
% Factors  : The factors in a cell array. Use any factors from 
%            any model. 
% G        : The core array. If 'G' is not defined it is assumed
%            that a PARAFAC model is being established.
%            Use G = [] in the PARAFAC case.
% Om       : Oblique mode.