pca_mean.m

一个利用Matlab实现PCA（主成分分析）

function [Y, lbd, cos2, CC] = pca_mean(X, P, labels)
%PCA Principal Components Analysis
%
%       Y = PCA(X, P) returns the first P principal components of X.
%       X must be a NxL matrix of N vectors in L dimensions.
%
%       Y = PCA(X, P, LABELS) displays several plots related to the PCA:
%         . projection of data on the first two principal axes,
%         . correlation circle,
%         . % inertia of the principal components,
%         . angles of the data w.r.t. the PC subspace.
%       If LABELS is a matrix containing N rows of text, each row is used
%       to label the corresponding data point, otherwise data go unlabelled.
%
%       Additional outputs: [Y, LBD, C2, CC] = pca_mean(X, P [, LABELS]), where
%       Y is the projection of X on the first P principal components,
%       LBD contains the % inertia of the required P components,
%       C2 is the squared cosine between projected data and original data,
%       CC is the correlation between the input features and the first P P.C.
%
% See also: .
if (nargin < 2) or (nargin > 4)
  error('PCANA: wrong number of arguments.')
end

[N, L] = size(X) ;
if (P < 1)
  error('PCANA: must return at least 1 principal component.')
end
if (P > L)
  error('PCANA: can''t get more principal components than dimensions.')
end

if (nargin == 2)
  verbose = 0 ;
end

%%%%%%%%%%

%R = corrcoef(X) ;
XM= ones(N,1) * mean(X) ;   % get the mean value
Xc = X - XM ;               % centalize the input data
[U, S, V] = svd(Xc, 0) ;    % singular decompose
%[V, S] = eig(R) ;
[l, idx] = sort(diag(S)') ; % order in descensivie order
l = fliplr(l) ;
idx = fliplr(idx) ;

% Coordonees centrees reduites.
XCR = (X - ones(N,1)*mean(X)) / diag(std(X)) ;  % be of no use, just neglect it
PC = (Xc * V(:,idx)' + XM) ;                    % get the principal component 
%PC = (XCR * V(:,idx) ) ;

% P first principal components:
Y = PC(:, 1:P) ;                               %  get the first P principal components

% P first portions of inertia:
lbd = l(1:P) / sum(l) ;                       % the contribution percentage of each component

% Squared cosines of projection vs. original data:
cos2 = sum((Y.^2)') ./ sum((XCR.^2)') ;           

% Correlations of features with PCs:
CC = V(:,idx) * diag(sqrt(l)) ;             % this valude will be dispayed as a curve , it indicates the egienvector of each components(notes that components
                                            % not principal components),
                                            % and  its value is
                                            % eigenvector*(sqrt(eigenvalue)
                                            % )
                                            

if (nargin == 3)
  figure
  subplot(2,2,1) ;
  plot(PC(:,1), PC(:,2), '+') ;
  if (size(labels, 1) == N)
    for i = 1:N
      text(PC(i,1), PC(i,2), [' ', labels(i,:)]) ;
    end
  end  
  xlabel('First principal axis') ;
  ylabel('Second principal axis') ;
  title('Data projection on the first two principal axes')
  subplot(2,2,2) ;
  plot(cos(2*pi*(0:100)/100), sin(2*pi*(0:100)/100)) ;
  hold on, plot([-1, 1], [0, 0]), plot([0, 0], [-1, 1]) ;
  plot(CC(:,1), CC(:,2), 'x') ;
  for k = 1:L
    text(CC(k,1)+.03, CC(k,2)+.03, num2str(k)) ;
  end
  xlabel('First PC') ;
  ylabel('Second PC') ;
  title('Correlation of features with first 2 PCs')
  subplot(2,2,3) ;
  hh = plot(l/sum(l), '--') ;
  set(hh, 'LineWidth', 3) ;
  hold on, hh = plot(cumsum(l)/sum(l)) ;
  set(hh, 'LineWidth', 3) ;
  axis([1 L 0 1]) ;
  plot([2 2], [0 1]) ;
  xlabel('PC no.') ;
  ylabel('Inertia') ;
  legend('PC inertia', 'Cumulated inertia') ;
  title('Repartition of inertia on the PCs') ;
  subplot(2,2,4) ;
  hold on
  hh = plot(acos(sqrt(cos2))/pi*180, 1:N, '+') ;
  set(hh, 'LineWidth', 3) ;
  hold on
  for i = 1:N
    hh = plot([0, acos(sqrt(cos2(i)))/pi*180], [i, i], '-') ;
    set(hh, 'LineWidth', 3) ;
  end
  set(gca, 'Box', 'on') ;
  set(gca, 'YLim', [1 N]) ;
  if (size(labels, 1) == N)
    set(gca, 'YTick', 1:N) ;
    set(gca, 'YTickLabels', labels) ;
  end
  xlabel('Angle (deg)')
  title('Angle of data with the PC plan')
end

CC = CC(:,1:P) ;