npls.m

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

   SSX=sum(sum(X.^2)); % To be used for evaluating the %var explained
end

% Check if weighting is tried used together with unimodality or orthogonality
if any(const==3)|any(const==1)
   if DoWeight==1
      disp(' ')
      disp(' Weighting is not possible together with unimodality and orthogonality.')
      disp(' It can be done using majorization, but has not been implemented here')
      disp(' Please contact the authors for further information')
      error
   end
end

% Acceleration
acc=-5;     
do_acc=1;   % Do acceleration every do_acc'th time
acc_pow=2;  % Extrapolate to the iteration^(1/acc_pow) ahead
acc_fail=0; % Indicate how many times acceleration have failed 
max_fail=4; % Increase acc_pow with one after max_fail failure
if showfit~=-1
   disp(' Line-search acceleration scheme initialized')
end

% Find initial guesses for the loadings if no initial values are given

% Use old loadings
if length(OldLoad)==ord % Use old values
   if showfit~=-1
      disp(' Using old values for initialization')
   end
   Factors=OldLoad;
   % Use DTLD
elseif Init==0
   if min(DimX)>1&ord==3&MissMeth==0
      if showfit~=-1
         disp(' Using direct trilinear decomposition for initialization')
      end
      [A,B,C]=dtld(reshape(X,DimX),Fac);
      A=real(A);B=real(B);C=real(C);
      % Check for signs and reflect if appropriate
      for f=1:Fac
         if sign(sum(A(:,f)))<0
            if sign(sum(B(:,f)))<0
               B(:,f)=-B(:,f);
               A(:,f)=-A(:,f);
            elseif sign(sum(C(:,f)))<0
               C(:,f)=-C(:,f);
               A(:,f)=-A(:,f);
            end
         end
         if sign(sum(B(:,f)))<0
            if sign(sum(C(:,f)))<0
               C(:,f)=-C(:,f);
               B(:,f)=-B(:,f);
            end
         end
      end
      Factors{1}=A;Factors{2}=B;Factors{3}=C;

   else
      if showfit~=-1
         disp(' Using singular values for initialization')
      end
      Factors=ini(X,Fac,2);
   end
   
   % Use SVD 
elseif Init==1
   if showfit~=-1
      disp(' Using singular values for initialization')
   end
   Factors=ini(X,Fac,2);
   
   % Use random (orthogonal)
elseif Init==2
   if showfit~=-1
      disp(' Using orthogonal random for initialization')
   end
   Factors=ini(X,Fac,1);
   
elseif Init==3
   error(' Initialization option set to three has been changed to 10')
   
   % Use several small ones of the above
elseif Init==10
   if showfit~=-1
      disp(' Using several small runs for initialization')
   end
   Opt=Options;
   Opt(5) = NaN;
   Opt(6) = NumbIteraInitia;
   Opt(2) = 0;
   ERR=[];
   [Factors,it,err] = parafac(X,Fac,Opt,const,[],[],Weights);
   ERR = [ERR;err];
   Opt(2) = 1;
   [F,it,Err] = parafac(X,Fac,Opt,const,[],[],Weights);
   ERR=[ERR;Err];
   if Err<err
      Factors=F;
      err=Err;
   end
   Opt(2)=2;
   for rep=1:3
      [F,it,Err]=parafac(X,Fac,Opt,const,[],[],Weights);
      ERR=[ERR;Err];
      if Err<err
         Factors=F;
         err=Err;
      end
   end
   if showfit~=-1
      disp(' ')
      disp(' Obtained fit-values')
      disp([' Method   Fit'])
      disp([' DTLD     ',num2str(ERR(1))])
      disp([' SVD      ',num2str(ERR(2))])
      disp([' RandOrth ',num2str(ERR(3))])
      disp([' RandOrth ',num2str(ERR(4))])
      disp([' RandOrth ',num2str(ERR(5))])
   end
else
   error(' Problem in PARAFAC initialization - Not set correct')
end
% Convert to old format
ff = [];
for f=1:length(Factors)
 ff=[ff;Factors{f}(:)];
end
Factors = ff;

% ALTERNATING LEAST SQUARES
err=SSX;
f=2*crit;
it=0;
connew=2;conold=1; % for missing values
ConstraintsNotRight = 0; % Just to ensure that iterations are not stopped if constraints are not yet fully imposed

if showfit~=-1
   disp(' ')
   disp(' Sum-of-Squares   Iterations  Explained')
   disp(' of residuals                 variation')
end

while ((f>crit) | (norm(connew-conold)/norm(conold)>MissConvCrit) | ConstraintsNotRight) & it<maxit
   conold=connew; % for missing values
   it=it+1;
   acc=acc+1; 
   if acc==do_acc;
      Load_o1=Factors;
   end
   if acc==do_acc+1;
      acc=0;Load_o2=Factors;
      Factors=Load_o1+(Load_o2-Load_o1)*(it^(1/acc_pow));
      % 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
         connew=model(id);
         errX=X-model;
         if DoWeight==0
            nerr=sum(sum(errX(idmiss2).^2));
         else
            nerr=sum(sum((Weights(idmiss2).*errX(idmiss2)).^2));
         end
      else
         if DoWeight==0
            nerr=sum(sum((X-model).^2));
         else
            nerr=sum(sum((X.*Weights-model.*Weights).^2));
         end
      end
      if nerr>err
         acc_fail=acc_fail+1;
         Factors=Load_o2;
         if acc_fail==max_fail,
            acc_pow=acc_pow+1+1;
            acc_fail=0;
            if showfit~=-1
               disp(' Reducing acceleration');
            end
         end
      else
         if MissMeth
            X(id)=model(id);
         end
      end
   end
   
   
   if DoWeight==0
      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=krb(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=krb(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=krb(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=krb(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