pca_hao.cpp

主成分分析程序

/*********************************/
/* Principal Components Analysis */
/*********************************/


#include <iostream>
#include <math.h>

using namespace std;

#define SIGN(a, b) ( (b) < 0 ? -fabs(a) : fabs(a) )

/**  Correlation matrix: creation  ***********************************/

void corcol(float** data, int n, int m, float** symmat)   

/* Create m * m correlation matrix from given n * m data matrix. */

{
float eps = 0.005;
float x, *mean, *stddev/*, *vector()*/;
int i, j, j1, j2;

/* Allocate storage for mean and std. dev. vectors */

mean = new float[m];
stddev = new float[m];

/* Determine mean of column vectors of input data matrix */

for (j = 0; j <= m-1; j++)
    {
    mean[j] = 0.0;
    for (i = 0; i <= n-1; i++)
        {
        mean[j] += data[i][j];
        }
    mean[j] /= (float)n;
    }

printf("\nMeans of column vectors:\n");
for (j = 0; j <= m-1; j++)  {
    printf("%7.2f",mean[j]);  }   printf("\n");

/* Determine standard deviations of column vectors of data matrix. */

for (j = 0; j <= m-1; j++)
    {
    stddev[j] = 0.0;
    for (i = 0; i <= n-1; i++)
        {
        stddev[j] += (   ( data[i][j] - mean[j] ) *
                         ( data[i][j] - mean[j] )  );
        }
        stddev[j] /= (float)n;
        stddev[j] = sqrt(stddev[j]);
        /* The following in an inelegant but usual way to handle
        near-zero std. dev. values, which below would cause a zero-
        divide. */
        if (stddev[j] <= eps) stddev[j] = 1.0;
    }

printf("\nStandard deviations of columns:\n");
for (j = 0; j <= m-1; j++) { printf("%7.1f", stddev[j]); }
printf("\n");

/* Center and reduce the column vectors. */

for (i = 0; i <= n-1; i++)
    {
    for (j = 0; j <= m-1; j++)
        {
        data[i][j] -= mean[j];
        x = sqrt((float)n);
        x *= stddev[j];
        data[i][j] /= x;
        }
    }

/* Calculate the m * m correlation matrix. */
for (j1 = 0; j1 <= m-1; j1++)
    {
    symmat[j1][j1] = 1.0;
    for (j2 = j1+1; j2 <= m-1; j2++)
        {
        symmat[j1][j2] = 0.0;
        for (i = 0; i <= n-1; i++)
            {
            symmat[j1][j2] += ( data[i][j1] * data[i][j2]);
            }
        symmat[j2][j1] = symmat[j1][j2];
        }
    }
    symmat[m-1][m-1] = 1.0;

	delete mean;
	delete stddev;


}

/**  Variance-covariance matrix: creation  *****************************/

void covcol(float** data, int n, int m, float** symmat)

/* Create m * m covariance matrix from given n * m data matrix. */

{
float *mean/*, *vector()*/;
int i, j, j1, j2;

/* Allocate storage for mean vector */

mean = new float[m];

/* Determine mean of column vectors of input data matrix */

for (j = 0; j <= m-1; j++)
    {
    mean[j] = 0.0;
    for (i = 0; i <= n-1; i++)
        {
        mean[j] += data[i][j];
        }
    mean[j] /= (float)n;
    }

printf("\nMeans of column vectors:\n");
for (j = 0; j <= m-1; j++)  {
    printf("%7.1f",mean[j]);  }   printf("\n");

/* Center the column vectors. */

for (i = 0; i <= n-1; i++)
    {
    for (j = 0; j <= m-1; j++)
        {
        data[i][j] -= mean[j];
        }
    }

/* Calculate the m * m covariance matrix. */
for (j1 = 0; j1 <= m-1; j1++)
    {
    for (j2 = j1; j2 <= m-1; j2++)
        {
        symmat[j1][j2] = 0.0;
        for (i = 0; i <= n-1; i++)
            {
            symmat[j1][j2] += data[i][j1] * data[i][j2];
            }
        symmat[j2][j1] = symmat[j1][j2];
        }
    }

delete mean;

}

/**  Sums-of-squares-and-cross-products matrix: creation  **************/

void scpcol(float** data, int n, int m, float**  symmat)

/* Create m * m sums-of-cross-products matrix from n * m data matrix. */

{
int i, j1, j2;

/* Calculate the m * m sums-of-squares-and-cross-products matrix. */

for (j1 = 0; j1 <= m-1; j1++)
    {
    for (j2 = j1; j2 <= m-1; j2++)
        {
        symmat[j1][j2] = 0.0;
        for (i = 0; i <= n-1; i++)
            {
            symmat[j1][j2] += data[i][j1] * data[i][j2];
            }
        symmat[j2][j1] = symmat[j1][j2];
        }
    }

}

/**  Reduce a real, symmetric matrix to a symmetric, tridiag. matrix. */

void tred2(float** a, int n, float* d,float* e)

/* Householder reduction of matrix a to tridiagonal form.
   Algorithm: Martin et al., Num. Math. 11, 181-195, 1968.
   Ref: Smith et al., Matrix Eigensystem Routines -- EISPACK Guide
        Springer-Verlag, 1976, pp. 489-494.
        W H Press et al., Numerical Recipes in C, Cambridge U P,
        1988, pp. 373-374.  */
{
int l, k, j, i;
float scale, hh, h, g, f;

for (i = n-1; i >= 1; i--)
    {
    l = i - 1;
    h = scale = 0.0;
    if (l > 0)
       {
       for (k = 0; k <= l; k++)
           scale += fabs(a[i][k]);
       if (scale == 0.0)
          e[i] = a[i][l];
       else
          {
          for (k = 0; k <= l; k++)
              {
              a[i][k] /= scale;
              h += a[i][k] * a[i][k];
              }
          f = a[i][l];
          g = f>0 ? -sqrt(h) : sqrt(h);
          e[i] = scale * g;
          h -= f * g;
          a[i][l] = f - g;
          f = 0.0;
          for (j = 0; j <= l; j++)
              {
              a[j][i] = a[i][j]/h;
              g = 0.0;
              for (k = 0; k <= j; k++)
                  g += a[j][k] * a[i][k];
              for (k = j+1; k <= l; k++)
                  g += a[k][j] * a[i][k];
              e[j] = g / h;
              f += e[j] * a[i][j];
              }
          hh = f / (h + h);
          for (j = 0; j <= l; j++)
              {
              f = a[i][j];
              e[j] = g = e[j] - hh * f;
              for (k = 0; k <= j; k++)
                  a[j][k] -= (f * e[k] + g * a[i][k]);
              }
         }
    }
    else
        e[i] = a[i][l];
    d[i] = h;
    }
d[0] = 0.0;
e[0] = 0.0;
for (i = 0; i <= n-1; i++)
    {
    l = i - 1;
	if(l>=0)
	{
    if (d[i])
       {
       for (j = 0; j <= l; j++)
           {
           g = 0.0;
           for (k = 0; k <= l; k++)
               g += a[i][k] * a[k][j];
           for (k = 0; k <= l; k++)
               a[k][j] -= g * a[k][i];
           }
       }
	}
       d[i] = a[i][i];
       a[i][i] = 1.0;
	if(l>=0)
       for (j = 0; j <= l; j++)
           a[j][i] = a[i][j] = 0.0;
	
	}
	
}

/**  Tridiagonal QL algorithm -- Implicit  **********************/

void tqli(float* d, float* e, int n,float** z)

{
int m, l, iter, i, k;
float s, r, p, g, f, dd, c, b;
//void erhand(char*);

for (i = 1; i <= n-1; i++)
    e[i-1] = e[i];
e[n-1] = 0.0;


for (l = 0; l <= n-1; l++)
    {
    iter = 0;
    do