jkparafac.m

Jackknife PARAFAC, 用在多维线性分解.

function [Factorsjk,Factorsini,ERR] = jkparafac(X,numfactors,Options,const,Weights,dibfac)

%JACK-KNIFE PARAFAC SEGMENTS OF DATA ARRAY X
% 
% The function has two modes: 1) normal jack-knifing and 2) a functionality for
% scaling and ordering sets of models. This functionality is built-in in 1) but is 
% useful to have if loadings have been generated in some other way.
%
% (1) function [Factorsjk,Factorsini,ERR] = jkparafac(X,numfactors,Options,const,Weights,dibfac)
% or
% (2) function [Fjkscaledandordered] = jkparafac(Factorsjk,Factorsini)
%
% ---------- 1: Do jack-knifing of PARAFAC model ----------
% Input
%
% X:            input array, which can be from three- to N-way
% numfactors:   number of factors of the PARAFAC model
% Options:      Optional parameters. If not given or set to zero or [], 
%               defaults will be used.
%                  Options(1) - Convergence criterion
%                  The relative change in fit for which the algorithm stops.
%                  Standard is 1e-6, but difficult data might require a
%                  lower value.
%                  Options(2) - Initialization method
%                     0  = fit using DTLD/GRAM for initialization
%                     1  = fit using SVD vectors for initialization
%                     2  = fit using random orthogonalized values for initialization
%                     10 = fit using the best-fitting models of several models
%                          fitted using a few iterations
%                  Options(3) - Plotting options (inside the parafac runs)
%                     0 = no plots
%                     1 = produces several graphical outputs
%                     2 = produces several graphical outputs (loadings also shown during iterations)
%                     3 = as 2 but no core consistency check (very slow for
%                         large arrays and/or many components) 
%                  Options(4) - Post-scaling of loadings
%                        0 or 1 = default scaling (columns in mode one carry the variance)
%                        2      = no scaling applied (hence fixed elements will not be modified
%                  Options(5) - How often to show fit
%                  Options(6) - Maximal number of iterations
%
% const:        A vector telling type of constraints put on the loadings of the
%               different modes. Same size as DimX but the i'th element tells
%               what constraint is on that mode.
%               0 => no constraint,
%               1 => orthogonality
%               2 => nonnegativity
%               3 => unimodality (and nonnegativitiy)
%               If const is not defined, no constraints are used.
%               For no constraints in a threeway problem const = [0 0 0]
%
% Weights       If a matrix of the same size as X is given, weighted regression
%               is performed using the weights in the matrix Weights. Statistically
%               the weights will usually contain the inverse error standard 
%               deviation of the particular element
%
% dibfac:       plotting options
%               1= plots of the jack-knife segments in the different modes
%               any other value= without plots
%
% Output
%
% Factorsjk: optimally scaled jack-knife segments of the PARAFAC loadings. For a
%            4 component solution to a 5 x 201 x 61 array the loadings Ajk, Bjk & Cjk
%            will be stored in a 3 element cell vector:
%            Factorsjk{1}=Ajk of size 5 x 4 x 5
%            Factorsjk{2}=Bjk of size 201 x 4 x 5
%            Factorsjk{3}=Cjk of size 61 x 4 x 5
%            etc.
%
%            The kth slab in the 3rd mode of Factorsjk{i} corresponds to the results
%            of the kth jack-knife segment for that mode. The kth row in the kth slab
%            of Factorsjk{1} is NaN because it corresponds to the segment in which
%            that object was left out
%           
%            Use FAC2LET.M for converting to "normal" output or simply extract the
%            components as e.g. Ajk = Factorsjk{1};
%
% Factorsini: loadings of the overall PARAFAC model. For a 4 component
%             solution to a 5 x 201 x 61 array the loadings Aini, Bini & Cini will be
%             stored in a 3 element cell vector:
%             Factorsini{1}=Aini of size 5 x 4
%             Factorsini{2}=Bini of size 201 x 4
%             Factorsini{3}=Cini of size 61 x 4
%             etc.
%             Use FAC2LET.M for converting to "normal" output or simply extract the
%             components as e.g. Aini = Factorsini{1};
%
% ERR: vector containing the fit of the model (the sum of squares of errors not including
%      missing elements) for the overall model and each jack-knife iteration
%
% ---------- 2: Scale and order sets of loaings to agree with overall model ----------
% Input
%
% Factorsjk: cell array containing loadings to be optimally scaled and permuted according
%            to Factorsini. For a 3-way array data the loadings Ajk, Bjk & Cjk should 
%            be stored in a 3-element cell vector:
%            Factorsjk{1}=Ajk of size i x number of factors x i
%            Factorsjk{2}=Bjk of size j x number of factors x i
%            Factorsjk{3}=Cjk of size k x number of factors x i
%            etc.
%            The kth row in the kth slab of Factorsjk{1} is NaN because it corresponds to
%            the segment in which that object was left out
%
% Factorsini: cell array containing the loadings of the overall PARAFAC model. For a 3-way
%             array data the loadings Aini, Bini & Cini should be stored in a 3-element
%             cell vector:
%             Factorsini{1}=Aini of size i x number of factors
%             Factorsini{2}=Bini of size j x number of factors
%             Factorsini{3}=Cini of size k x number of factors
%
% Output
%
% Fjkscaledandordered: optimally scaled and ordered segments of Factorsjk
%
% Please refer to this m-file through
%
%     Jordi Riu and Rasmus Bro, Jack-knife technique for outlier detection
%     and estimation of standard errors in PARAFAC models, Chemometrics and
%     Intelligent Laboratory Systems 65 (2003) 35-49
%
%  I/0:
%     (1) function [Factorsjk,Factorsini,ERR] = jkparafac(X,numfactors,Options,Weights,dibfac)
%  or
%     (2) function [Fjkscaledandordered] = jkparafac(Factorsjk,Factorsini)
%
% see ripplot, impplot


% 2004, march, made jkparafac selfcontained so that it runs without the
% nway toolbox
% 2004, april, replaced fminunc with fminsearch so it runs without the
% optimization toolbox
% 2004, july, fixed minor bug that caused occasional crashes in some
% situations
% 2005, Sept, minor update due to updates in nway toolbox
% 2006, Feb, minor bug in plotting fixed (Thanks to Carina Rubingh)

jkp=0;

global dim
global Ag;
global Ajk;
global numfac

% check if we only have to perform scale and permutation (jkp=1)

[res,lletres1]=size(class(X));
cc1=class(X);

[res,lletres2]=size(class(numfactors));
cc2=class(numfactors);

if lletres1==4 & cc1(1)=='c' & cc1(2)=='e' & lletres2==4 & cc2(1)=='c' & cc2(2)=='e'
    jkp=1;
    Factorsjkper=X;
    Factorsiniper=numfactors;
    [res,dim]=size(Factorsiniper);
    for i=1:dim;
        [DimX(i) numfac]=size(Factorsiniper{i});
    end
    ERR=[];
    Options=[];
    dibfac=0;
end;

if exist('dibfac')~=1;
    dibfac=0;
end;

if exist('Options')~=1;
    Options=[];
end;

if exist('const')~=1;
    const=[];
end;

if exist('Weights')~=1;
    Weights=[];
end;


if jkp==0;
    DimX=size(X);
    [res,dim]=size(size(X));

    ERR=[];
    Factorsjk=cell(1,dim);
    Aunf=cell(1,1);

    numfac=numfactors;

% Global model

    disp(' ')    
    disp ('          Global model')
    disp(' ')
    [Factorsini,it,err]=parafac(X,numfac,Options,const,[],[],Weights);
    ERR=[ERR err];
end;
% Jack-knife segments / scaling and permutation problems

jks=1;

for jkit=1:jks:DimX(1);
    if jkp==0;
    
        disp(' ')    
        disp (['          segment number ' num2str(jkit)])
        disp(' ')
    end;
    
        if jkit==1
            if jkp==0;
                xnew=reshape(X(jks+1:DimX(1),:),[(DimX(1)-jks) DimX(2:end)]);
                if ~isempty(Weights)&length(size(Weights))==length(DimX)
                  Weightsnew=reshape(Weights(jks+1:DimX(1),:),[(DimX(1)-jks) DimX(2:end)]);
                else
                  Weightsnew = [];
                end
                aininew=Factorsini{1}(jks+1:DimX(1),:);
            end;
            segment=[jks+1:DimX(1)];
        elseif jkit>1
            if jkp==0;
                x1=reshape(X(1:jkit-1,:),[(jkit-1) DimX(2:end)]);
                x2=reshape(X(jkit+jks:DimX(1),:),[(DimX(1)-jkit-jks+1) DimX(2:end)]);
                
                aininew1=Factorsini{1}(1:jkit-1,:);
                aininew2=Factorsini{1}(jkit+jks:DimX(1),:);
                xnew=[x1; x2];
                aininew=[aininew1; aininew2];
                if ~isempty(Weights)&length(size(Weights))==length(DimX)
                  w1=reshape(Weights(1:jkit-1,:),[(jkit-1) DimX(2:end)]);
                  w2=reshape(Weights(jkit+jks:DimX(1),:),[(DimX(1)-jkit-jks+1) DimX(2:end)]);
                  Weightsnew=[w1;w2];
                else
                  Weightsnew = [];
                end
                
            end;
            segment=[1:jkit-1 jkit+jks:DimX(1)];
        end;
    
    if jkp==0;
        initialfactors={aininew Factorsini{2:end}};
        [factorsnew,it,err]=parafac(xnew,numfac,Options,const,initialfactors,[],Weightsnew);
        ERR=[ERR err];
    end;

    
% solving the scaling and permutation problems

    for i=1:dim;
        if jkp==1;
            Factorsini{i}=Factorsiniper{i};
            if i==1;
                factorsnew{i}=Factorsjkper{i}(segment,:,jkit);
            elseif i~=1
                factorsnew{i}=Factorsjkper{i}(:,:,jkit);
            end;
        end;
        if i==1;
            [reca,perma]=RecFac('MitBur',factorsnew{i},Factorsini{i}(segment,:));
            [factorsnew{i}]=FacPerm(perma,factorsnew{i});
            Tg{i}=Factorsini{i}(segment,:);
            Tjk{i}=factorsnew{i};
        elseif i~=1
            [reca,perma]=RecFac('MitBur',factorsnew{i},Factorsini{i});
            [factorsnew{i}]=FacPerm(perma,factorsnew{i});
            Tg{i}=Factorsini{i};
            Tjk{i}=factorsnew{i};    
        end;
        Op=optimset('Display','off','LargeScale','off');
    end;
    Ag=Tg;
    Ajk=Tjk;  
    
    
    %[ST,valT]=fminunc(@MinTjkpar,ones(numfac*dim,1),Op);    
    [ST,valT]=fminsearch(@MinTjkpar,ones(numfac*dim,1),Op);
    
    for i=1:dim;
        res=factorsnew{i}*diag(ST(numfac*(i-1)+1:i*numfac));
        if i==1 & jkit==1
            res=[NaN*ones(1,numfac);res];
        elseif i==1 & jkit>1 & jkit<DimX(1)
            res=[res(1:jkit-1,:);NaN*ones(1,numfac);res(jkit:DimX(1)-jks,:)];
        elseif i==1 & jkit==DimX(1)
            res=[res;NaN*ones(1,numfac)];            
        end;
        Factorsjk{i}(:,:,jkit)=res;
    end;       
end;

% number of rows and columns in each subplot

rows=1;
columns=1;
while rows*columns<dim
    columns=columns+1;
    if rows*columns<dim
    rows=rows+1;
    end;
end;

if dibfac==1
    for j=1:numfac
        figure
        ara=squeeze(Factorsjk{1}(:,j,:));
        subplot(rows,columns,1),
        plot(ara,'o')
        title(['mode number 1, component number ' num2str(j)])
        
        for k=2:dim;
            ara=squeeze(Factorsjk{k}(:,j,:));
            subplot (rows,columns,k),plot(ara)
            title(['mode number ' num2str(k) ', component number ' num2str(j)])
        end;
    end;
end;

function [Recovery,Perm] = RecFac(Meas,varargin)
if isodd(nargin-1)
   error
end
Dims = (nargin - 1) / 2;
switch Meas
case {'MaxCor'}
   d = ppp(varargin{Dims},varargin{Dims-1});
   e = ppp(varargin{Dims*2},varargin{Dims * 2 - 1});
   for i = Dims - 2:-1:1
      d = ppp(d,varargin{i});
      e = ppp(e,varargin{Dims + i});
   end
   for o = 1:size(e,2)
      Co = 0;
      for p = 1:size(d,2)
         CCL = corrcoef([d(:,p),e(:,o)]);
         if abs(CCL(2)) > Co;
            Co = abs(CCL(2));
            Pos  = p;
         end
      end
      Perm(Pos,o) = 1;
      Recovery(o) = Co;
   end
case 'MaxCos'
   for o = 1:size(varargin{Dims+1},2)
      Co = 0;
      for p = 1:size(varargin{1},2)
         CCL = 1;
         for q = 1:Dims
            CCL = CCL .* (varargin{q}(:,p)' * varargin{Dims + q}(:,o))/(norm(varargin{q}(:,p))*norm(varargin{Dims + q}(:,o)));
         end
         if CCL > Co;
            Co = CCL;
            Pos  = p;
         end
      end
      Perm(Pos,o) = 1;
      Recovery(o) = Co;
   end
case 'MitBur'
    
   if size(varargin{1},2) < size(varargin{Dims + 1},2)
      error('The model with the largest n. of factors must come first')
   end
   Factors1 = varargin(1:Dims);
   Factors2 = varargin(Dims + 1:2 * Dims);   
   Rk1      = size(Factors1{1},2);
   Rk2      = size(Factors2{1},2);
   R        = perms(1:Rk1);
   Perm     = zeros(Rk1);
   Sig      = ones(1,size(varargin{Dims + 1},2));
   Per      = [];
   Co       = -inf;
   for o = 1:size(R,1)
      CCL = ones(Rk2);
      for q = 1:Dims
         Cross          = normit(Factors1{q}(:,R(o,1:Rk2)))' * normit(Factors2{q});
         temp_sign(:,q) = sign(diag(Cross));
         CCL            = CCL .* Cross;
      end
      Rec = mean(diag(CCL));
      if Rec > Co
         Co       = Rec;
         Recovery = diag(CCL);
         Sig      = temp_sign;
         Per      = o;
      end
   end   
   if ~isempty(Per)
      Perm(sub2ind([Rk1 Rk1],R(Per,:),1:Rk1)) = 1;
   else 
      Perm(1:Rk1,1:Rk1) = eye(Rk1);
      Recovery          = NaN*zeros(1,Rk1);
   end
   Recovery = Recovery(:);
   
case 'MSE'
   Rk   = size(varargin{Dims + 1},2);
   Recovery = NaN * ones(Dims,Rk);
   for p = 1:Rk
      M = 1;
      X = 1;
      for q = Dims:-1:1
         Recovery(q,p) = 100 * (1 - sum((varargin{q}(:,p) - varargin{Dims + q}(:,p)).^2) / sum(varargin{Dims + q}(:,p).^2));
      end
   end
   Perm = eye(Rk);
end
if ~strcmp(Meas,'MSE')
   Recovery = Recovery(:);
end

function varargout = FacPerm(Perm,varargin)

for i=1:nargin-1
   varargout{i} = varargin{i} * Perm;
end

function b = isodd(a)
if rem(a,2)
   b = logical(1);
else