C*********************** Contents ****************************************
C* Principal Components Analysis: C, 638 lines. ****************************
C* Sample input data set (final 36 lines). *********************************
C***************************************************************************
/*********************************/
/* Principal Components Analysis */
/*********************************/
/*********************************************************************/
/* Principal Components Analysis or the Karhunen-Loeve expansion is a
classical method for dimensionality reduction or exploratory data
analysis. One reference among many is: F. Murtagh and A. Heck,
Multivariate Data Analysis, Kluwer Academic, Dordrecht, 1987.
Author:
F. Murtagh
Phone: + 49 89 32006298 (work)
+ 49 89 965307 (home)
Earn/Bitnet: fionn@dgaeso51, fim@dgaipp1s, murtagh@stsci
Span: esomc1::fionn
Internet: [email protected]
F. Murtagh, Munich, 6 June 1989 */
/*********************************************************************/
#include <stdio.h>
#include <string.h>
#include <math.h>
#define SIGN(a, b) ( (b) < 0 ? -fabs(a) : fabs(a) )
main(argc, argv)
int argc;
char *argv[];
{
FILE *stream;
int n, m, i, j, k, k2;
float **data, **matrix(), **symmat, **symmat2, *vector(), *evals, *interm;
void free_matrix(), free_vector(), corcol(), covcol(), scpcol();
void tred2(), tqli();
float in_value;
char option, *strncpy();
/*********************************************************************
Get from command line:
input data file name, #rows, #cols, option.
Open input file: fopen opens the file whose name is stored in the
pointer argv[argc-1]; if unsuccessful, error message is printed to
stderr.
*********************************************************************/
if (argc != 5)
{
printf("Syntax help: PCA filename #rows #cols option\n\n");
printf("(filename -- give full path name,\n");
printf(" #rows \n");
printf(" #cols -- integer values,\n");
printf(" option -- R (recommended) for correlation analysis,\n");
printf(" V for variance/covariance analysis\n");
printf(" S for SSCP analysis.)\n");
exit(1);
}
n = atoi(argv[2]); /* # rows */
m = atoi(argv[3]); /* # columns */
strncpy(&option,argv[4],1); /* Analysis option */
printf("No. of rows: %d, no. of columns: %d.\n",n,m);
printf("Input file: %s.\n",argv[1]);
if ((stream = fopen(argv[1],"r")) == NULL)
{
fprintf(stderr, "Program %s : cannot open file %s\n",
argv[0], argv[1]);
fprintf(stderr, "Exiting to system.");
exit(1);
/* Note: in versions of DOS prior to 3.0, argv[0] contains the
string "C". */
}
/* Now read in data. */
data = matrix(n, m); /* Storage allocation for input data */
for (i = 1; i <= n; i++)
{
for (j = 1; j <= m; j++)
{
fscanf(stream, "%f", &in_value);
data[i][j] = in_value;
}
}
/* Check on (part of) input data.
for (i = 1; i <= 18; i++) {for (j = 1; j <= 8; j++) {
printf("%7.1f", data[i][j]); } printf("\n"); }
*/
symmat = matrix(m, m); /* Allocation of correlation (etc.) matrix */
/* Look at analysis option; branch in accordance with this. */
switch(option)
{
case 'R':
case 'r':
printf("Analysis of correlations chosen.\n");
corcol(data, n, m, symmat);
/* Output correlation matrix.
for (i = 1; i <= m; i++) {
for (j = 1; j <= 8; j++) {
printf("%7.4f", symmat[i][j]); }
printf("\n"); }
*/
break;
case 'V':
case 'v':
printf("Analysis of variances-covariances chosen.\n");
covcol(data, n, m, symmat);
/* Output variance-covariance matrix.
for (i = 1; i <= m; i++) {
for (j = 1; j <= 8; j++) {
printf("%7.1f", symmat[i][j]); }
printf("\n"); }
*/
break;
case 'S':
case 's':
printf("Analysis of sums-of-squares-cross-products");
printf(" matrix chosen.\n");
scpcol(data, n, m, symmat);
/* Output SSCP matrix.
for (i = 1; i <= m; i++) {
for (j = 1; j <= 8; j++) {
printf("%7.1f", symmat[i][j]); }
printf("\n"); }
*/
break;
default:
printf("Option: %s\n",option);
printf("For option, please type R, V, or S\n");
printf("(upper or lower case).\n");
printf("Exiting to system.\n");
exit(1);
break;
}
/*********************************************************************
Eigen-reduction
**********************************************************************/
/* Allocate storage for dummy and new vectors. */
evals = vector(m); /* Storage alloc. for vector of eigenvalues */
interm = vector(m); /* Storage alloc. for 'intermediate' vector */
symmat2 = matrix(m, m); /* Duplicate of correlation (etc.) matrix */
for (i = 1; i <= m; i++) {
for (j = 1; j <= m; j++) {
symmat2[i][j] = symmat[i][j]; /* Needed below for col. projections */
}
}
tred2(symmat, m, evals, interm); /* Triangular decomposition */
tqli(evals, interm, m, symmat); /* Reduction of sym. trid. matrix */
/* evals now contains the eigenvalues,
columns of symmat now contain the associated eigenvectors. */
printf("\nEigenvalues:\n");
for (j = m; j >= 1; j--) {
printf("%18.5f\n", evals[j]); }
printf("\n(Eigenvalues should be strictly positive; limited\n");
printf("precision machine arithmetic may affect this.\n");
printf("Eigenvalues are often expressed as cumulative\n");
printf("percentages, representing the 'percentage variance\n");
printf("explained' by the associated axis or principal component.)\n");
printf("\nEigenvectors:\n");
printf("(First three; their definition in terms of original vbes.)\n");
for (j = 1; j <= m; j++) {
for (i = 1; i <= 3; i++) {
printf("%12.4f", symmat[j][m-i+1]); }
printf("\n"); }
/* Form projections of row-points on first three prin. components. */
/* Store in 'data', overwriting original data. */
for (i = 1; i <= n; i++) {
for (j = 1; j <= m; j++) {
interm[j] = data[i][j]; } /* data[i][j] will be overwritten */
for (k = 1; k <= 3; k++) {
data[i][k] = 0.0;
for (k2 = 1; k2 <= m; k2++) {
data[i][k] += interm[k2] * symmat[k2][m-k+1]; }
}
}
printf("\nProjections of row-points on first 3 prin. comps.:\n");
for (i = 1; i <= n; i++) {
for (j = 1; j <= 3; j++) {
printf("%12.4f", data[i][j]); }
printf("\n"); }
/* Form projections of col.-points on first three prin. components. */
/* Store in 'symmat2', overwriting what was stored in this. */
for (j = 1; j <= m; j++) {
for (k = 1; k <= m; k++) {
interm[k] = symmat2[j][k]; } /*symmat2[j][k] will be overwritten*/
for (i = 1; i <= 3; i++) {
symmat2[j][i] = 0.0;
for (k2 = 1; k2 <= m; k2++) {
symm
