ncrossdecompn.m

多维数据分析,有nPLS,PARAFAC,TURKER等

      
      
      % Xvalidated Model of data
      ModelXval = zeros(I,J)*NaN;
      for s = 1:Segments
         disp([' Segment ',num2str(s),'/',num2str(Segments),' - Comp. ',num2str(possibleCombs(f1,1:end-1)),'/',num2str(FacMax)])
         %GT
         Xnow    = permute(zeros(DimX),Perm);
         dimsadd = DimX;
         for j = 1:length(DimX)
             dimsadd(j) = delta(j);
             Xnow       = cat(j,Xnow,zeros(dimsadd));
             index2{j}  = 1:DimX(j);
             dimsadd(j) = dimsadd(j) + DimX(j);
         end
         Xnow(s:Segments:end)     = NaN;
         Xnow                     = reshape(permute(Xnow(index2{:}),PermI),DimX(1),prod(DimX(2:end)));
         Pos                      = isnan(Xnow);
         [M,Mean]                 = decomp(Method,Xnow + X,DimX,possibleCombs(f1,1:end-1),s,Segments,Cent,I,J,Param);
         model                    = M + ones(I,1)*Mean';
         ModelXval(Pos)           = model(Pos);
         %GT end
      end
      
      XvalResult.Xval(f1,:) = [100*(1 - sum( (X(id) - ModelXval(id)).^2)/sum(OffsetCorrectedData(id).^2)) possibleCombs(f1,1:end-1)];
      XvalResult.XvalModel{f1} =  ModelXval;
   end
   
end



if Show&FacMin-FacMax~=0
   if Method(1:3) == 'pca'
      Nam = 'PCA';
   elseif Method(1:3) == 'tuc'
      Nam = 'Tucker';
   elseif Method(1:3) == 'par'
      Nam = 'PARAFAC';
   elseif Method(1:3) == 'nip'
      Nam = 'NIPALS';
   end
   
   figure      
   save jjj
   if lower(Method(1:3))~='tuc'
      bar(FacMin:FacMax,[XvalResult.Fit(FacMin:FacMax,1) XvalResult.Xval(FacMin:FacMax,1)],.76,'grouped')
   else
      % extract the ones with lowest Xval fit (for each # total comp) for plotting
      fx = [];
      f5 =[];
      for f1 = 1:max(possibleCombs(:,end))
         f2 = find(possibleCombs(:,end)==f1);
         if length(f2)
            [f3,f4] = max(XvalResult.Xval(f2));
            f5 = [f5,f2(f4)];
         end
      end
      fx = [possibleCombs(f5,end) XvalResult.Fit(f5,1) XvalResult.Xval(f5,1)];
      bar(fx(:,1),fx(:,2:3),.76,'grouped')
      for f1 = 1:size(fx,1)
         f6=text(fx(f1,1),95,['[',num2str(possibleCombs(f5(f1),1:end-1)),']']);
         set(f6,'Rotation',270)
      end
   end
   
   g=get(gca,'YLim');
   set(gca,'YLim',[max(-20,g(1)) 100])
   legend('Fitted','Xvalidated',0)
   titl = ['Xvalidation results (',Nam,')'];
   if Cent
      titl = [titl ,' - centering'];
   else
      titl = [titl ,' - no centering'];
   end
   title(titl,'FontWeight','Bold')
   xlabel('Total number of components')
   ylabel('Percent variance explained')
end   



function [M,Mean,parameters] = decomp(Method,X,DimX,f,s,Segments,Cent,I,J,parameters);

Conv = 0;
it = 0;
maxit = 500;
% Initialize
if Cent
   Mean = nanmean(X)';
else
   Mean = zeros(J,1);
end

if lower(Method(1:3)) == 'par'
   Xc = reshape(X- ones(I,1)*Mean',DimX);
   if exist('parameters')==1
      fact = parafac(Xc,f,[1e-5 10 0 0 NaN maxit],[],parameters.fact);
   else
      fact = parafac(Xc,f,[1e-5 10 0 0 NaN maxit]);
   end
   M = reshape(nmodel(fact),DimX(1),prod(DimX(2:end)));
   parameters.fact=fact;
elseif lower(Method(1:3)) == 'tuc'
   Xc = reshape(X- ones(I,1)*Mean',DimX);
   if exist('parameters')==1
      [fact,G] = tucker(Xc,f,[1e-2 0 0 0 NaN maxit],[],[],parameters.fact,parameters.G);
   else
      [fact,G] = tucker(Xc,f,[1e-2 0 0 0 NaN maxit]);
   end
   parameters.fact=fact;
   parameters.G=G;
   M = reshape(nmodel(fact,G),DimX(1),prod(DimX(2:end)))      ;
elseif lower(Method) == 'pca'|lower(Method) == 'nip'
   Xc = reshape(X- ones(I,1)*Mean',DimX(1),prod(DimX(2:end)));
   [t,p] = pcanipals(X- ones(I,1)*Mean',f,0);
   parameters.t=t;
   parameters.p=p;
   M = t*p';
else
   error(' Name of method not recognized') 
end
Fit = X - M - ones(I,1)*Mean';
Fit = sum(Fit(find(~isnan(X))).^2);

% Iterate
while ~Conv
   it     = it+1;
   FitOld = Fit;
   
   % Fit multilinear part
   Xcent = X - ones(I,1)*Mean';
   
   if Method(1:3) == 'par'
      fact = parafac(reshape(Xcent,DimX),f,[1e-2 0 0 0 NaN maxit],[],fact);
      M = reshape(nmodel(fact),DimX(1),prod(DimX(2:end)));
   elseif Method(1:3) == 'tuc'
      [fact,G] = tucker(reshape(Xcent,DimX),f,[1e-2 0 0 0 NaN maxit],[0 0 0],zeros(size(G)),fact,G);
      M = reshape(nmodel(fact,G),DimX(1),prod(DimX(2:end)));
   elseif Method == 'pca'
      [t,p] = pcals(Xcent,f,0,t,p,0);
      M = t*p';
   elseif Method == 'nip'
      [t,p] = pcanipals(Xcent,f,0);
      M = t*p';
   end
   
   % Find offsets
   if Cent
      x = X;
      mm=M+ones(I,1)*Mean';
      x(find(isnan(X)))=mm(find(isnan(X)));
      Mean = mean(x)';
   end
   
   
   %Find fit
   Fit = X - M - ones(I,1)*Mean';
   Fit = sum(Fit(find(~isnan(X))).^2);
   if abs(Fit-FitOld)/FitOld<1e-8 | it > 1500
      Conv = 1;
   end
   
end
disp([' Fit ',num2str(Fit),' using ',num2str(it),' it.'])



function [t,p] = pcals(X,F,cent,t,p,show);

%  LEAST SQUARES PCA WITH MISSING ELEMENTS
%  20-6-1999
% 
%  Calculates a least squares PCA model. Missing elements 
%  are denoted NaN. The solution is NOT nested, so one has
%  to calculate a new model for each number of components.


ShowMeFitEvery = 20;
MaxIterations  = 5;
[I,J]=size(X);
Xorig      = X;
Miss       = find(isnan(X));
NotMiss    = find(~isnan(X));
m          = t*p';
X(Miss)    = m(Miss);
ssX    = sum(X(NotMiss).^2);

Fit    = 3;
OldFit = 6;
it     = 0;

while abs(Fit-OldFit)/OldFit>1e-3 & it < MaxIterations;
   it      = it +1;
   OldFit  = Fit;
   
   [t,s,p] = svds(X,F);
   t       = t*s;
   
   Model   = t*p';
   X(Miss) = Model(Miss);
   Fit     = sum(sum( (Xorig(NotMiss) - Model(NotMiss)).^2));
   
   if ~rem(it,ShowMeFitEvery)&show
      disp(['    Fit after ',num2str(it),' it. :',num2str(RelFit),'%'])
   end
   
end


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

% NIPALS-PCA WITH MISSING ELEMENTS
% cent: One if centering is to be included, else zero


[I,J]=size(X);
rand('state',sum(100*clock))

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      = rand(I,1);
   P      = rand(J,1);
   Fit    = 2;
   FitOld = 3;
   
   while abs(Fit-FitOld)/FitOld>1e-7 & it < 100;
      FitOld  = Fit;
      it      = it +1;
      
      for j = 1:J
         id=find(NotMiss(:,j));
         if length(id)==0
            id,end
         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

function Xc = nanmean(X)

if isempty(X)
   Xc = NaN;
   return
end

i = isnan(X);
j = find(i);
i = sum(i);
X(j) = 0;
Num = size(X,1)-i;
Xc = sum(X);
i = find(Num);
Xc(i) = Xc(i)./Num(i);
Xc(find(~Num))=NaN;