1 files changed, 79 insertions, 0 deletions
diff --git a/modules/statistics/macros/princomp.sci b/modules/statistics/macros/princomp.sci
new file mode 100755
index 000000000..b7f3e1cf8
--- /dev/null
+++ b/modules/statistics/macros/princomp.sci
@@ -0,0 +1,79 @@
+// Scilab ( http://www.scilab.org/ ) - This file is part of Scilab
+// Copyright (C) 2000 - INRIA - Carlos Klimann
+//
+// This file must be used under the terms of the CeCILL.
+// This source file is licensed as described in the file COPYING, which
+// you should have received as part of this distribution.  The terms
+// are also available at
+// http://www.cecill.info/licences/Licence_CeCILL_V2.1-en.txt
+//
+
+function [facpr,comprinc,lambda,tsquare]=princomp(x,eco)
+    //
+    //This  function performs  several  computations known  as
+    //"principal component  analysis".
+    //The  idea  behind this  method  is  to  represent in  an
+    //approximative  manner a  cluster of  n individuals  in a
+    //smaller  dimensional subspace.  In order  to do  that it
+    //projects the cluster onto a subspace.  The choice of the
+    //k-dimensional projection subspace is  made in such a way
+    //that  the distances  in  the projection  have a  minimal
+    //deformation: we are looking for a k-dimensional subspace
+    //such that the squares of the distances in the projection
+    //is  as  big  as  possible  (in  fact  in  a  projection,
+    //distances can only stretch).  In other words, inertia of
+    //the projection  onto the k dimensional  subspace must be
+    //maximal.
+    //
+
+    //References: Saporta, Gilbert, Probabilites,  Analyse des
+    //Donnees et Statistique, Editions Technip, Paris, 1990.
+    //
+    //derived from pca,
+    //  author: carlos klimann
+    //
+
+    [lhs, rhs] = argn(0);
+    if rhs == 0 then
+        error(msprintf(gettext("%s: Wrong number of input arguments: %d or %d expected.\n"),"princomp",1,2))
+    end
+    if rhs<2 then eco=%f,end
+    if type(eco)<>4 then
+        error(msprintf(_("%s: Wrong type for input argument #%d: Boolean expected."),"princomp",2))
+    end
+
+    if x==[] then
+        facpr=[];
+        comprinc=[];
+        lambda=[];
+        tsquare=[];
+        return;
+    end
+
+    [rowx,colx]=size(x)
+
+    //center data
+    y=center(x,1)
+
+    //compute eigenvectors of  y'*y using svd
+    if eco then
+        [U,lambda,facpr]=svd(y,"e");
+    else
+        [U,lambda,facpr]=svd(y);
+    end
+    lambda=diag(lambda).^2/(rowx-1);//normalize
+    comprinc=y*facpr
+    if rowx <= colx then
+        if eco then
+            lambda(rowx)=[]
+            facpr(:,rowx)=[];
+            comprinc(:,rowx) = [];
+        else
+            lambda(rowx:$)=[]
+            comprinc(:,rowx:$) = [];
+        end
+    end
+    q=find(lambda<=max(rowx,colx)*%eps*lambda(1),1)
+    if q==[] then q=size(lambda,"*"),end
+    tsquare=(rowx-1)*sum(U(:,1:q).^2,2)
+endfunction