parafac.m

强大的多维工具箱.应用在Matlab中,可分析多纬数据结构.直接安装.

    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/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


% 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

% 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(:);

  %   % Instead of above, do signs so as to make them as "natural" as possible
  %   Factors = signswtch(Factors,reshape(X,DimX));
  %   DIDN't WORK (TOOK AGES FOR 7WAY DATA)


  if showfit~=-1
    disp(' Components have been reflected according to convention')
  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;


if Plt==1|Plt==2|Plt==3
  %   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

function swloads = signswtch(loads,X);

%SIGNSWTCH switches sign of multilinear models so that signs are in
%accordance with majority of data
%
%
% I/O swloads = signswtch(loads,X);
%
% Factors must be a cell with the loadings. If Tucker or NPLS, then the
% last element of the cell must be the core array

try % Does not work in older versions of matlab
  warning('off','MATLAB:divideByZero');
end
sizeX=size(X);
order = length(sizeX);
for i=1:order;
  F(i) = size(loads{i},2);
end


if isa(X,'dataset')% Then it's a SDO
  inc=X.includ;
  X = X.data(inc{:});
end


% Compare centered X with center loading vector

if length(loads)==order % PARAFAC
  % go through each component and then update in the end
  for m = 1:order % For each mode determine the right sign
    for f=1:F(1) % one factor at the time
      s=[];
      a = loads{m}(:,f);
      x = permute(X,[m 1:m-1 m+1:order]);
      for i=1:size(x(:,:),2); % For each column
        id = find(~isnan(x(:,i)));
        if length(id)>1
          try
            c = corrcoef(x(id,i),a(id));
          catch
            disp('Oops - something wrong in signswtch - please send a note to [e-mail]')
            whos
          end
          if isnan(c(2,1))
            s(i)=0;
          else
            s(i) = c(2,1)*length(id); % Weigh correlation by number of elements so many-miss columns don't influence too much
          end
        else
          s(i) = 0;
        end
      end
      S(m,f) = sum(s);
    end
  end

  % Use S to switch signs. If the signs of S (for each f) multiply to a
  % positive number the switches are performed. If not, the mode of the
  % negative one with the smallest absolute value is not switched.

  for f = 1:F(1)
    if sign(prod(S(:,f)))<1 % Problem: make the smallest negative positive to avoid switch of that
      id = find(S(:,f)<0);
      [a,b]=min(abs(S(id,f)));
      S(id(b(1)),f)=-S(id(b(1)),f);
    end
  end
  % Now ok, so switch what needs to be switched
  for f = 1:F(1)
    for m = 1:order
      if sign(S(m,f))<1
        loads{m}(:,f)=-loads{m}(:,f);
      end
    end
  end




elseif length(loads)==(order+1) % NPLS/Tucker

  % go through each mode and update and correct core accordinglu
  for m = 1:order % For each mode determine the right sign
    for f=1:F(m) % one factor at the time
      a = loads{m}(:,f);
      x = permute(X,[m 1:m-1 m+1:order]);
      for i=1:size(x(:,:),2); % For each column
        id = find(~isnan(x(:,i)));
        if length(id)>1
          c = corrcoef(x(id,i),a(id));
          if isnan(c(2,1))
            s(i)=0;
          else
            s(i) = c(2,1)*length(id); % Weigh correlation by number of elements so many-miss columns don't influence too much
          end
        else
          s(i) = 0;
        end
      end
      if sum(s) < 0
        % turn around
        loads{m}(:,f) = -loads{m}(:,f);

        % Then switch the core accordingly
        G = loads{order+1};
        G = permute(G,[m 1:m-1 m+1:order]);
        sizeG = size(G);
        G = reshape(G,sizeG(1),prod(sizeG)/sizeG(1));
        G(f,:) = -G(f,:);
        G = reshape(G,sizeG);
        G = ipermute(G,[m 1:m-1 m+1:order]);
        loads{order+1} = G;
      end
    end
  end


else
  error('Unknown model type in SIGNS.M')
end

swloads = loads;