npls.m

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

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(krb(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*krb(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


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

% NIPALS-PCA WITH MISSING ELEMENTS
% 20-6-1999
%
% Calculates a NIPALS PCA model. Missing elements 
% are denoted NaN. The solution is nested
% 
% Comparison for data with missing elements
% NIPALS : Nested    , not least squares, not orthogonal solutoin
% LSPCA  : Non nested, least squares    , orthogonal solution
% 
% I/O
% [t,p,Mean,Fit,RelFit] = pcanipals(X,F,cent);
% 
% X   : Data with missing elements set to NaN
% F   : Number of componets
% cent: One if centering is to be included, else zero
% 
% Copyright
% Rasmus Bro
% KVL 1999
% [e-mail]
%

[I,J]=size(X);
if any(sum(isnan(X))==I)|any(sum(isnan(X)')==J)
   error(' One column or row only contains missing')
end

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      = nanmean(X')';
   P      = nanmean(X)';
   Fit    = 2;
   FitOld = 3;
   
   while abs(Fit-FitOld)/FitOld>1e-7 & it < 1000;
      FitOld  = Fit;
      it      = it +1;
      
      for j = 1:J
         id=find(NotMiss(:,j));
         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

Model   = t*p' + ones(I,1)*Mean;
Fit     = sum(sum( (Xorig(find(NotMiss)) - Model(find(NotMiss))).^2));
RelFit  = 100*(1-Fit/ssX);

function y = nanmean(x)
if isempty(x) % Check for empty input.
    y = NaN;
    return
end
nans = isnan(x);
i = find(nans);
x(i) = zeros(size(i));

if min(size(x))==1,
  count = length(x)-sum(nans);
else
  count = size(x,1)-sum(nans);
end
i = find(count==0);
count(i) = ones(size(i));
y = sum(x)./count;
y(i) = i + NaN;