%
%  PCA - Principal Components Analysis
%
%        [ t, p, pc ] = pca( x, dpc );
%  where
%        x   =   data set (rows are points and
%                columns are variables)
%        dpc =   desired percentage of variance
%                to be explained
%        t   =   columns are score vectors
%        p   =   columns are loading vectors
%        pc  =   percentage of variance explained
%                by each component
%
%  Calculation of the principal component score and loading
%  vectors for a data set by the NIPALS algorithm
%
%  S.Wold, "Cross validatory estimation of the number of
%   components in factor and principal components models",
%   Technometrics, 20, 397-406, (1978)
%
%  Bob Newell, November 1994
%
function [ t, p, pc ] = pca( x, percent )
[nr nc] = size( x );
%
% scale and centre column-wise
%
ord = [];
for i = 1:nc
  m = mean( x(:,i) );
  s = std( x(:,i) );
  ord = [ ord s ];
  x(:,i) = ( x(:,i) - m ) / s;
end
varo = norm( x ); varl = varo;
na = 0; conv = 0; t = []; p = []; pc = [];
%
% start of the loop
%
while ( conv == 0 )
%
% component number
%
ord
%  pause
  na = na + 1;
%
% find column with largest variance
%
  [s i] = max( ord );
  tv = x(:,i);
%
% normalise the chosen column
%
  tv = tv / norm( tv );
  d = 1;
  n = 0;
  while ( d > 1e-10 )
%
%   calculate p and then a new t
%
    pv = ( tv' * x / ( tv' * tv ) )';
    pn = pv / sqrt( pv' * pv );
    tn = x * pn / ( pn' * pn );
%    [ tv  tn ]
%
%   inner convergence check 
%
    n = n + 1;
    if n > 25
      d = 1e-20;
    else
      d = ( tn - tv )' * ( tn - tv ) / nr / ( tn' * tn );
    end
%    pause
    tv = tn;
  end
  t = [ t tv ];
  p = [ p pn ];
%
% calculate residuals for next iteration
%
  x = x - tv * pn';
  var = norm( x );
  pc = [ pc ( 100 * ( varl - var ) / varo ) ];
  varl = var;
  ord(i) = 0;
%
% outer convergence check
%
  if ( na == nc ) 
    conv = 1;
  end
  if 100 * ( varo - var ) / varo > percent
    conv = 1;
  end
%
% end of main loop
%
end
