octave_packages/statistics-1.1.3/princomp.m

   1 ## Author: Paul Kienzle <pkienzle@users.sf.net>
   2 ## This program is granted to the public domain.
   3
   4 ## -*- texinfo -*-
   5 ## @deftypefn {Function File} {[@var{pc}, @var{z}, @var{w}, @var{Tsq}] =} princomp (@var{X})
   6 ##
   7 ## Compute principal components of @var{X}.
   8 ##
   9 ## The first output argument @var{pc} is the principal components of @var{X}.
  10 ## The second @var{z} is the transformed data, and @var{w} is the eigenvalues of
  11 ## the covariance matrix of @var{X}. @var{Tsq} is the Hotelling's @math{T^2}
  12 ## statistic for the transformed data.
  13 ## @end deftypefn
  14
  15 function [pc,z,w,Tsq] = princomp(X)
  16   C = cov(X);
  17   [U,D,pc] = svd(C,1);
  18   if nargout>1, z = center(X)*pc; end
  19   if nargout>2, w = diag(D); end
  20   if nargout>3, Tsq = sumsq(zscore(z),2);
  21     warning('XXX FIXME XXX Tsq return from princomp fails some tests');
  22   end
  23 endfunction
  24 %!shared pc,z,w,Tsq,m,x
  25
  26 %!test
  27 %! x=[1,2,3;2,1,3]';
  28 %! [pc,z,w,Tsq]=princomp(x);
  29 %! m=[sqrt(2),sqrt(2);sqrt(2),-sqrt(2);-2*sqrt(2),0]/2;
  30 %! m(:,1) = m(:,1)*sign(pc(1,1));
  31 %! m(:,2) = m(:,2)*sign(pc(1,2));
  32
  33 %!assert(pc,m(1:2,:),10*eps);
  34 %!assert(z,-m,10*eps);
  35 %!assert(w,[1.5;.5],10*eps);
  36 %!assert(Tsq,[4;4;4]/3,10*eps);
  37
  38 %!test
  39 %! x=x';
  40 %! [pc,z,w,Tsq]=princomp(x);
  41 %! m=[sqrt(2),sqrt(2),0;-sqrt(2),sqrt(2),0;0,0,2]/2;
  42 %! m(:,1) = m(:,1)*sign(pc(1,1));
  43 %! m(:,2) = m(:,2)*sign(pc(1,2));
  44 %! m(:,3) = m(:,3)*sign(pc(3,3));
  45
  46 %!assert(pc,m,10*eps);
  47 %!assert(z(:,1),-m(1:2,1),10*eps);
  48 %!assert(z(:,2:3),zeros(2),10*eps);
  49 %!assert(w,[1;0;0],10*eps);
  50 %!xtest
  51 %! assert(Tsq,1,10*eps);