Daniel@0: function [P,V,me,l] = pcaproj(D,arg1,arg2)
Daniel@0: 
Daniel@0: %PCAPROJ Projects data vectors using Principal Component Analysis.
Daniel@0: %
Daniel@0: % [P,V,me,l] = pcaproj(D, odim)
Daniel@0: % P =          pcaproj(D, V, me)
Daniel@0: %
Daniel@0: %  Input and output arguments ([]'s are optional)
Daniel@0: %   D      (matrix) size dlen x dim, the data matrix
Daniel@0: %          (struct) data or map struct            
Daniel@0: %   odim   (scalar) how many principal vectors are used
Daniel@0: %  
Daniel@0: %   P      (matrix) size dlen x odim, the projections
Daniel@0: %   V      (matrix) size dim x odim, principal eigenvectors (unit length)
Daniel@0: %   me     (vector) size 1 x dim, center point of D
Daniel@0: %   l      (vector) size 1 x odim, the corresponding eigenvalues, 
Daniel@0: %                   relative to total sum of eigenvalues
Daniel@0: %                   
Daniel@0: % See also SAMMON, CCA.
Daniel@0: 
Daniel@0: % Contributed to SOM Toolbox 2.0, February 2nd, 2000 by Juha Vesanto
Daniel@0: % Copyright (c) by Juha Vesanto
Daniel@0: % http://www.cis.hut.fi/projects/somtoolbox/
Daniel@0: 
Daniel@0: % juuso 191297 070200
Daniel@0: 
Daniel@0: %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
Daniel@0: 
Daniel@0: error(nargchk(2, 3, nargin)); % check the number of input arguments
Daniel@0: 
Daniel@0: % the data
Daniel@0: if isstruct(D), 
Daniel@0:   if strcmp(D.type,'som_map'), D=D.codebook; else D=D.data; end
Daniel@0: end
Daniel@0: [dlen dim] = size(D);
Daniel@0: 
Daniel@0: if nargin==2, 
Daniel@0: 
Daniel@0:   odim = arg1;
Daniel@0:     
Daniel@0:   % autocorrelation matrix
Daniel@0:   A = zeros(dim);
Daniel@0:   me = zeros(1,dim);
Daniel@0:   for i=1:dim, 
Daniel@0:     me(i) = mean(D(isfinite(D(:,i)),i)); 
Daniel@0:     D(:,i) = D(:,i) - me(i); 
Daniel@0:   end  
Daniel@0:   for i=1:dim, 
Daniel@0:     for j=i:dim, 
Daniel@0:       c = D(:,i).*D(:,j); c = c(isfinite(c));
Daniel@0:       A(i,j) = sum(c)/length(c); A(j,i) = A(i,j); 
Daniel@0:     end
Daniel@0:   end
Daniel@0:   
Daniel@0:   % eigenvectors, sort them according to eigenvalues, and normalize
Daniel@0:   [V,S]   = eig(A);
Daniel@0:   eigval  = diag(S);
Daniel@0:   [y,ind] = sort(abs(eigval)); 
Daniel@0:   eigval  = eigval(flipud(ind));
Daniel@0:   V       = V(:,flipud(ind)); 
Daniel@0:   for i=1:odim, V(:,i) = (V(:,i) / norm(V(:,i))); end
Daniel@0:   
Daniel@0:   % take only odim first eigenvectors
Daniel@0:   V = V(:,1:odim);
Daniel@0:   l = abs(eigval)/sum(abs(eigval));
Daniel@0:   l = l(1:odim); 
Daniel@0: 
Daniel@0: else % nargin==3, 
Daniel@0: 
Daniel@0:   V = arg1;
Daniel@0:   me = arg2;
Daniel@0:   odim = size(V,2);    
Daniel@0:   D = D-me(ones(dlen,1),:);
Daniel@0:   
Daniel@0: end
Daniel@0:   
Daniel@0: % project the data using odim first eigenvectors
Daniel@0: P = D*V;
Daniel@0: 
Daniel@0: %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%