/* PRINCIPAL COMPONENTS
Submitted 12 April 1996 by Mico Loretan (loretanm@frb.gov)
The algorithm is based on discussion in
Theil's _Principles of Econometrics_.
The code is provided free for public non-commercial use.
There are no performance guarantees.

Mico Loretan
loretanm@frb.gov

Remark: Principal Components Analysis
can be interpreted as a special case of
factor analysis, where all the "factors"
are mutually uncorrelated and
constructed to sequentially maximize
"explanatory power" (in the R^2 sense.)
*/

/* PROGRAM CODE: PRINCOMP.SRC */
/*
:: princomp.src
::
:: Mico Loretan (loretanm@frb.gov), February 1996
::
:: For a matix X of order (nxk), calculate
::  - The first "j" principal components (matrix of order nxj), in descending
::    order of the amount of variance "explained" (j<=k)
::  - Proportion of variance explained by each of the j components
::  - "Factor loadings" of columns of x  with respect to each
::    of the principal components.
:: Note: This algorithm is based on: Henri Theil, Priciples of
:: Econometrics, 1971, New York: Wiley, pp. 46-56 (principal comps.)
*/

proc (3)=princomp(x,j);
    /*
    :: Inputs:  x  data matrix, of order nxk (assume: n>k, full rank)
    ::          j  number of principal components to
    ::             be computed (j<=k)
    :: Outputs: p  matrix of the first j principal components of X (nxj)
    ::          v  vector of fractions of variance explained (jx1)
    ::          a  matrix of factor loadings, such that x = p*a+error
    ::             (matrix of order jxk)
    :: Notes: (i) This proc makes use of Gauss's proc "eigrs2()," which
    :: computes the eigenvalues _and_ the normalized eigenvectors of a
    :: real symmetric matrix, in ascending order of the eigenvalues.
    :: (ii) The principal components are computed indirectly, by
    :: first computing the eigenvectors of the x'x matrix.  Why?
    :: Computing the xx' matrix (and its eigenvectors) directly
    :: will exceed available memory for moderate and large values of n!
    */


    local k, first_j, eigvals, eigvecs, evals, evecs,
        p, v, a;               /* three output variables */

    k = cols(x);

    /* quick input sanity check */
    if j>k;
        errorlog "Can't compute more than cols(x) principal components.";
        end;
    endif;


    /* indices of first j principal components */
    first_j = seqa(k,-1,j);


    /* First, compute eigenvectors of x'x */

    {eigvals,eigvecs} = eigrs2(x'x);    /* sorted in ascending order */

    evals = submat(eigvals,first_j,1);  /* pick off j largest values */
    evecs = submat(eigvecs,0,first_j);  /* pick off corr. eigenvectors */


    /*
    :: We need the eigenvectors of x'x to have non-unit length; in
    :: fact, the length of each eigenvector must equal the square root
    :: of the corresponding eigenvalue.
    :: Note: (a') is also equal to the matrix of "factor loadings"
    :: for the j principal components (Theil, p. 54).
    :: (I.e., we have x=p*a'+v, where (a') is the matrix of
    :: least squares coefficients of the regression of x on p.
    :: If j=k, then x=p*a', and all variation in x is "explained.")
    */

    a = evecs .* sqrt(evals');


    /*
    :: Second, get eigenvectors of xx' (= principal components of x!)
    */

    p = (x*a)./(evals');    /* these eigenvectors have unit length */


    /*
    :: Third, compute fraction of total variation explained by
    :: each of the j principal components.
    */

    v = evals/sumc(eigvals);


retp(p,v,a');
endp; /* princomp() */