ncrossdecomp.m

多维数据处理：MATLAB源程序用于处理多维数据

function XvalResult = ncrossdecomp(Method,X,DimX,FacMin,FacMax,Segments,Cent,Show);

% CROSSVALIDATION FOR MULTI-WAY DECOMPOSITION MODELS
%
% $ Version 1.0301 $ Date 28. June 1999 $ Not compiled $
%
% See also:
% 'ncross' 
%
% Copyright, 1999 - 
% 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]
%
% This file performs cross-validation of decomposition models
% PARAFAC, PCA, and Tucker. The cross-validation is performed
% such that part of the data are set to missing, the model is 
% fitted to the remaining data, and the residuals between fitted
% and true left-out elements is calculated. This is performed
% 'Segments' times such that all elements are left out once.
% The segments are chosen by taking every 'Segments' element of
% X(:), i.e. from the vectorized array. If X is of size 5 x 7, 
% and three segemnts are chosen ('Segments' = 3), then in the 
% first of three models, the model is fitted to the matrix
% 
% |x 0 0 x 0 0 x|
% |0 x 0 0 x 0 0|
% |0 0 x 0 0 x 0|
% |x 0 0 x 0 0 x|
% |0 x 0 0 x 0 0|
% 
% where x's indicate missing elements. After fitting the residuals
% in the locations of missing values are calculated. After fitting
% all three models, all residuals have been calculated.
% 
% Note that the number of segments must be chosen such that no columns
% or rows contain only missing elements (the algorithm will check this).
% Using 'Segments' = 7, 9, or 13 will usually achieve that.
% 
% I/O
% XvalResult = ncrossdecomp(Method,X,DimX,FacMin,FacMax,Segments,Cent,Show);
% 
% INPUT
% Method   : 'parafac', 'tucker', 'pca', or 'nipals'
%            For Tucker only Tucker3 models with equal
%            number of components is currently available.
%            For PCA the least squares model is calculated.
%            Thus, offsets and parameters are calculated in
%            a least squares sense unlike the method NIPALS,
%            which calculates the PCA model using an ad hoc 
%            approach for handling missing data (as in 
%            standard chemometric software).
% X        : Data array unfolded (to an I x JK matrix in case of
%            three-way data)
% DimX     : Dimension of X ( = [I J K] in case of three-way)
% FacMin   : 
% FacMax   : 
% Segments : 
% Cent     : 
% Show     : 
%
% OUTPUT
% XvalResult 
%
%	Copyright
%	Rasmus Bro 1997
%	Denmark
%	E-mail [e-mail]


% uses NANSUM,

%
if 4==2
   Fmin = 1;
   Fmax = 5;
   load brod;
   for i=1:10
      respca= ncrossdecomp('pca',X(:,[1:2 4:31 33:end]),DimX,Fmin,Fmax,7,1,1);save anal;
      eval(['respca',num2str(i),'=respca;']);
   end
   for i=1:10
      respf= ncrossdecomp('par',X,DimX,Fmin,Fmax,7,1,1);save anal;
      eval(['respf',num2str(i),'=respf;']);
   end
   for i=1:10
      restuc= ncrossdecomp('tuc',X,DimX,Fmin,Fmax,7,1,1);save anal;
      eval(['restuc',num2str(i),'=restuc;']);
   end
   for i=1:10
      resnip= ncrossdecomp('nip',X(:,[1:2 4:31 33:end]),DimX,Fmin,Fmax,7,1,1);save anal;
      eval(['resnip',num2str(i),'=resnip;']);
   end

end

[I,J] = size(X);

if exist('Show')~=1
   Show = 1;
end


% Check if the selected segmentation works (does not produce rows/columns of only missing)
out = ones(I,J);
out(1:Segments:end)=NaN;
out(find(isnan(X))) = NaN;
if any(sum(isnan(out))==I)
   error(' The chosen segmentation leads to columns of only missing elements')
elseif any(sum(isnan(out'))==J)
   error(' The chosen segmentation leads to rows of only missing elements')
end
   

XvalResult.Fit = (1:FacMax)'*NaN;
XvalResult.Xval = (1:FacMax)'*NaN;
for f = 1:FacMin-1
   XvalResult.XvalModel{f} =  'Not fitted';
   XvalResult.FittedModel{f} = 'Not fitted';
end

for f = FacMin:FacMax
   
   % Fitted model
   disp([' Total model - Comp. ',num2str(f),'/',num2str(FacMax)])
   [M,Mean,Param] = decomp(Method,X,DimX,f,1,Segments,Cent,I,J);
   Model = M + ones(I,1)*Mean';
   id    = find(~isnan(X));
   OffsetCorrectedData = X - ones(I,1)*Mean';
   XvalResult.Fit(f) = 100*(1 - sum( (X(id) - Model(id)).^2)/sum(OffsetCorrectedData(id).^2));
   XvalResult.FittedModel{f} = Model;

   
   % Xvalidated Model of data
   ModelXval = zeros(I,J)*NaN;
   for s = 1:Segments
      disp([' Segment ',num2str(s),'/',num2str(Segments),' - Comp. ',num2str(f),'/',num2str(FacMax)])
      Xnow = X;
      Xnow(s:Segments:end) = NaN;
      [M,Mean] = decomp(Method,Xnow,DimX,f,s,Segments,Cent,I,J,Param);
      model = M + ones(I,1)*Mean';
      ModelXval(s:Segments:end) = model(s:Segments:end);
   end
   
   XvalResult.Xval(f) = 100*(1 - sum( (X(id) - ModelXval(id)).^2)/sum(OffsetCorrectedData(id).^2));
   XvalResult.XvalModel{f} =  ModelXval;
end

if Show&FacMin-FacMax~=0
   if Method(1:3) == 'pca'
      Nam = 'PCA';
   elseif Method(1:3) == 'tuc'
      Nam = 'Tucker';
   elseif Method(1:3) == 'par'
      Nam = 'PARAFAC';
   elseif Method(1:3) == 'nip'
      Nam = 'NIPALS';
   end
   
   figure      
   bar(FacMin:FacMax,[XvalResult.Fit(FacMin:FacMax) XvalResult.Xval(FacMin:FacMax)],.76,'grouped')
   g=get(gca,'YLim');
   set(gca,'YLim',[max(-20,g(1)) 100])
   legend('Fitted','Xvalidated',0)
   title(['Xvalidation results (',Nam,')'],'FontWeight','Bold')
   xlabel('Number of components')
   ylabel('Percent variance explained')
end   



function [M,Mean,parameters] = decomp(Method,X,DimX,f,s,Segments,Cent,I,J,parameters);

   Conv = 0;
   it = 0;
   
   % Initialize
   if Cent
      Mean = nanmean(X)';
   else
      Mean = zeros(J,1);
   end
   
   if Method(1:3) == 'par'
      if exist('parameters')==1
         fact = parafac(X- ones(I,1)*Mean',DimX,f,[1e-5 10 0 0 NaN 20],[],parameters.fact);
      else
         fact = parafac(X- ones(I,1)*Mean',DimX,f,[1e-5 10 0 0 NaN 20]);
      end
      M = nmodel(fact,DimX,f);
      parameters.fact=fact;
   elseif Method(1:3) == 'tuc'
      if exist('parameters')==1
         [fact,G] = tucker(X- ones(I,1)*Mean',DimX,[f f f],[1e-2 0 0 0 NaN 20],[],[],parameters.fact,parameters.G);
      else
         [fact,G] = tucker(X- ones(I,1)*Mean',DimX,[f f f],[1e-2 0 0 0 NaN 20]);
      end
      parameters.fact=fact;
      parameters.G=G;
      M = nmodel(fact,G,DimX,[f f f]);
   elseif Method == 'pca'|Method == 'nip'
      [t,p] = pcanipals(X- ones(I,1)*Mean',f,0);
      parameters.t=t;
      parameters.p=p;
      M = t*p';
   end
   Fit = X - M - ones(I,1)*Mean';
   Fit = sum(Fit(find(~isnan(X))).^2);
      
   % Iterate
   while ~Conv
      it     = it+1;
      FitOld = Fit;
      
      % Fit multilinear part
      Xcent = X - ones(I,1)*Mean';

      if Method(1:3) == 'par'
         fact = parafac(Xcent,DimX,f,[1e-2 0 0 0 NaN 20],[],fact);
         M = nmodel(fact,DimX,f);
      elseif Method(1:3) == 'tuc'
         [fact,G] = tucker(Xcent,DimX,[f f f],[1e-2 0 0 0 NaN 20],[0 0 0],zeros(size(G)),fact,G);
         M = nmodel(fact,G,DimX,[f f f]);
      elseif Method == 'pca'
         [t,p] = pcals(Xcent,f,0,t,p,0);
         M = t*p';
      elseif Method == 'nip'
         [t,p] = pcanipals(Xcent,f,0);
         M = t*p';
      end
      
      % Find offsets
      if Cent
         x = X;
         mm=M+ones(I,1)*Mean';
         x(find(isnan(X)))=mm(find(isnan(X)));
         Mean = mean(x)';
      end
      
      
      %Find fit
      Fit = X - M - ones(I,1)*Mean';
      Fit = sum(Fit(find(~isnan(X))).^2);
      if abs(Fit-FitOld)/FitOld<1e-8 | it > 1500
         Conv = 1;
      end

   end
   disp([' Fit ',num2str(Fit),' using ',num2str(it),' it.'])
   

   
function [t,p] = pcals(X,F,cent,t,p,show);
   
%  LEAST SQUARES PCA WITH MISSING ELEMENTS
%  20-6-1999
% 
%  Calculates a least squares PCA model. Missing elements 
%  are denoted NaN. The solution is NOT nested, so one has
%  to calculate a new model for each number of components.


ShowMeFitEvery = 20;
MaxIterations  = 5;
[I,J]=size(X);
Xorig      = X;
Miss       = find(isnan(X));
NotMiss    = find(~isnan(X));
m          = t*p';
X(Miss)    = m(Miss);
ssX    = sum(X(NotMiss).^2);

Fit    = 3;
OldFit = 6;
it     = 0;

while abs(Fit-OldFit)/OldFit>1e-3 & it < MaxIterations;
   it      = it +1;
   OldFit  = Fit;
 
   [t,s,p] = svds(X,F);
   t       = t*s;
   
   Model   = t*p';
   X(Miss) = Model(Miss);
   Fit     = sum(sum( (Xorig(NotMiss) - Model(NotMiss)).^2));
   
   if ~rem(it,ShowMeFitEvery)&show
      disp(['    Fit after ',num2str(it),' it. :',num2str(RelFit),'%'])
   end

end
 
 
function [t,p,Mean] = pcanipals(X,F,cent);

% NIPALS-PCA WITH MISSING ELEMENTS
% cent: One if centering is to be included, else zero
 

[I,J]=size(X);
rand('state',sum(100*clock))

Xorig      = X;
Miss       = isnan(X);
NotMiss    = ~isnan(X);
ssX    = sum(X(find(NotMiss)).^2);
Mean   = zeros(1,J);
if cent
   Mean    = nanmean(X);
end
X      = X - ones(I,1)*Mean;

t=[];
p=[];

for f=1:F
   Fit    = 3;
   OldFit = 6;
   it     = 0;
   T      = rand(I,1);
   P      = rand(J,1);
   Fit    = 2;
   FitOld = 3;
   
   while abs(Fit-FitOld)/FitOld>1e-7 & it < 100;
      FitOld  = Fit;
      it      = it +1;
      
      for j = 1:J
         id=find(NotMiss(:,j));
if length(id)==0
            id,end
         P(j) = T(id)'*X(id,j)/(T(id)'*T(id));
      end
      P = P/norm(P);
      
      for i = 1:I
         id=find(NotMiss(i,:));
         T(i) = P(id)'*X(i,id)'/(P(id)'*P(id));
      end
      
      Fit = X-T*P';
      Fit = sum(Fit(find(NotMiss)).^2);
   end
   t = [t T];
   p = [p P];
   X = X - T*P';

end

function Xc = nanmean(X)

if isempty(X)
   Xc = NaN;
   return
end

i = isnan(X);
j = find(i);
i = sum(i);
X(j) = 0;
Num = size(X,1)-i;
Xc = sum(X);
i = find(Num);
Xc(i) = Xc(i)./Num(i);
Xc(find(~Num))=NaN;