corcond.m

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

function [Consistency,G,stdG,Target]=corcond(X,DimX,Fac,Factors,Weights,Plot);
% $ Version 1.02 $ Date 28. July 1998 $ Not compiled $
%
% See also:
% 'unimodal' 'monreg' 'fastnnls'
%
% 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]
%
%
% CORe CONsistency DIAgnostics (corcondia)
% Performs cocote of a PARAFAC model and returns the cocote plot
% as well as the degree of consistency (100 % is max).
%
% Consistency=corcond(X,DimX,Fac,Factors,Weights);
% 
% INPUT
% X        : Data array unfolded to I x JK matrix
% DimX     : Dimensions of array [I J K]
% Fac      : Number of factors
% Factors  : Factors given in standard format [A(:);B(:);C(:)];
% Weights  : Optional weights (otherwise skip input or give an empty array [])
% Plot     = 0 or not given => no plots are produced
%          = 1              => normal corcondia plot
%          = 2              => corcondia plot with standard deviations 
%
% OUTPUT
% The core consistency given as the percentage of variation in a Tucker3 core
% array consistent with the theoretical superidentity array. Max value is 100%
% Consistencies well below 70-90% indicates that either too many components
% are used or the model is otherwise mis-specified.
% 
%	Copyright
%	Rasmus Bro 1997
%	Denmark
%	E-mail [e-mail]

if nargin<6
  Plot=0;
end
if nargin<5
  Weights=0;
end

ord=length(DimX);
l_idx=0;
for i=1:ord
  l_idx=[l_idx sum(DimX(1:i))*Fac];
end

if sum(DimX*Fac)~=length(Factors)
  error(' Size of Factors does not fit the dimension and number of factors given')
end


% Scale all loadings to same magnitude
magn=ones(Fac,1);
for i=1:ord
   L=reshape(Factors(l_idx(i)+1:l_idx(i+1)),DimX(i),Fac);
   for f=1:Fac
     magn(f)=magn(f)*norm(L(:,f));
     L(:,f)=L(:,f)/norm(L(:,f));
   end
   Factors(l_idx(i)+1:l_idx(i+1))=L(:);
end
% Magn holds the singular value of each component. Scale each loading vector by 
% the cubic root (if three-way) so all loadings of a component have the same variance

magn = magn.^(1/ord);
for i=1:ord
   L=reshape(Factors(l_idx(i)+1:l_idx(i+1)),DimX(i),Fac);
   for f=1:Fac
     L(:,f)=L(:,f)*magn(f);
   end
   Factors(l_idx(i)+1:l_idx(i+1))=L(:);
end


% Make diagonal array holding the magnitudes
Ident=nident(Fac,ord);
DimIdent=ones(1,ord)*Fac;
Ident=nshape(Ident,DimIdent,ord);

% Make matrix of Kronecker product of all loadings expect the large; Z = kron(C,B ... )
  NewFac=[];
  NewFacNo=[];
  for i=ord:-1:1
    Z=reshape(Factors(l_idx(i)+1:l_idx(i+1)),DimX(i),DimIdent(i));
    % Check its of full rank or adjust core and use less columns
    rankZ=rank(Z);
    if rankZ<Fac
       %OLD out=Z(:,rankZ+1:Fac);Z=Z(:,1:rankZ);H=[[eye(rankZ)] pinv(Z)*out];Ident=H*Ident;
       [q,r]=qr(Z);
       Ident=r*Ident;
       Z=q;
       DimIdent(i)=size(r,1);
    end
    if i>1
      Ident=nshape(Ident,DimIdent([i:ord 1:i-1]),ord);
    end
    NewFac=[Z(:);NewFac];
    NewFacNo=[rankZ NewFacNo];
  end
if length(NewFac)~=length(Factors)
  disp(' ')
  disp(' Some loading matrices do not have full column rank. Therefore')
  disp(' the model has been restated in performing CORCONDIA. This does')
  disp(' not affect the model')
end
Factors=NewFac;
Fac=NewFacNo;

[G,stdG]=regg(X,DimX,Fac,Factors,Weights);

 Ident=Ident(:);
 Target=Ident;
 [a,b]=sort(abs(Ident));
 b=flipud(b);
 Ident=Ident(b);
 GG=G(b);
 stdGG=stdG(b);
 bNonZero=find(Ident);
 bZero=find(~Ident);

 ssG=sum(G(:).^2);
 Consistency=100*(1-sum((Target-G).^2)/ssG);

 if Plot
    clf
    Ver=version;
    Ver=Ver(1);
    if Fac>1
       eval(['set(gcf,''Name'',''Diagonality test'');']);
       if Ver>4
          plot([Ident(bNonZero);Ident(bZero)],'y','LineWidth',3)
          hold on
          plot(GG(bNonZero),'ro','LineWidth',3)
          plot(length(bNonZero)+1:prod(Fac),GG(bZero),'gx','LineWidth',3)
          if Plot==2
            line([[1:length(G)];[1:length(G)]],[GG GG+stdGG]','LineWidth',1,'Color',[0 0 0])
            line([[1:length(G)];[1:length(G)]],[GG GG-stdGG]','LineWidth',1,'Color',[0 0 0])
          end
          hold off
          title(['Core consistency ',num2str(Consistency),'% (yellow target)'],'FontWeight','bold','FontSize',12)          
       else
          plot([Ident(bNonZero);Ident(bZero)],'y')
          hold on
          plot(GG(bNonZero),'ro')
          plot(length(bNonZero)+1:prod(Fac),GG(bZero),'gx')
          if Plot==2
            line([[1:length(G)];[1:length(G)]],[GG GG+stdGG]','LineWidth',1,'Color',[0 0 1])
            line([[1:length(G)];[1:length(G)]],[GG GG-stdGG]','LineWidth',1,'Color',[0 0 1])
          end
          hold off
          title(['Core consistency ',num2str(Consistency),'% (yellow target)'])
       end
       xlabel('Core elements (green should be zero/red non-zero)')
       ylabel('Core Size')
    else
       eval(['set(gcf,''Name'',''Diagonality test'');']);
       title(['Core consistency ',num2str(Consistency),'% (yellow target)'])
       xlabel('Core elements (green should be zero/red non-zero)')
       ylabel('Size')
       plot(GG(bNonZero),'ro')
       title(['Core consistency ',num2str(Consistency),'%'])
       xlabel('Core elements (red non-zero)')
       ylabel('Core Size')
    end
 end