pca_hao.cpp

主成分分析程序

      {
      for (m = l; m <= n-2; m++)
          {
          dd = fabs(d[m]) + fabs(d[m+1]);
          if (fabs(e[m]) + dd == dd) break;
          }
          if (m != l)
             {
             g = (d[l+1] - d[l]) / (2.0 * e[l]);
             r = sqrt((g * g) + 1.0);
             g = d[m] - d[l] + e[l] / (g + SIGN(r, g));
             s = c = 1.0;
             p = 0.0;
             for (i = m-1; i >= l; i--)
                 {
                 f = s * e[i];
                 b = c * e[i];
                 if (fabs(f) >= fabs(g))
                    {
                    c = g / f;
                    r = sqrt((c * c) + 1.0);
                    e[i+1] = f * r;
                    c *= (s = 1.0/r);
                    }
                 else
                    {
                    s = f / g;
                    r = sqrt((s * s) + 1.0);
                    e[i+1] = g * r;
                    s *= (c = 1.0/r);
                    }
                 g = d[i+1] - p;
                 r = (d[i] - g) * s + 2.0 * c * b;
                 p = s * r;
                 d[i+1] = g + p;
                 g = c * r - b;
                 for (k = 0; k <= n-1; k++)
                     {
                     f = z[k][i+1];
                     z[k][i+1] = s * z[k][i] + c * f;
                     z[k][i] = c * z[k][i] - s * f;
                     }
                 }
                 d[l] = d[l] - p;
                 e[l] = g;
                 e[m] = 0.0;
             }
          }  while (m != l);
      }
 }



main()
{
int  n, m,  i, j, k, k2;
float **data, /***matrix(),*/ **symmat, **symmat2, /**vector(), */*evals, *interm, t;
//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.
*********************************************************************/


    cin>>m;
    cin>>n;
	data = new float*[n];
	for(i=0;i<=n-1;i++)
		data[i] = new float[m];/* Storage allocation for input data */

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

    /* 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 = new float*[m];
	for(i=0;i<=m-1;i++)
		symmat[i] = new float[m];/* Allocation of correlation (etc.) matrix */

   /* Look at analysis option; branch in accordance with this. */

	cin>>option;

     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 = new float[m];     /* Storage alloc. for vector of eigenvalues */
    interm = new float[m];    /* Storage alloc. for 'intermediate' vector */
    //symmat2 = new[] float*(m);
	symmat2=new float*[m];
	for(i=0;i<=m-1;i++)
		symmat2[i] = new float[m];/* Duplicate of correlation (etc.) matrix */
    for (i = 0; i <= m-1; i++) 
	{
     for (j = 0; j <= m-1; j++) 
	 {
      symmat2[i][j] = symmat[i][j]; /* Needed below for col. projections */
                             
	 }
                            
	}

    tred2(symmat, m, evals, interm);  /* Triangular decomposition */
    for(i=0;i<=m-1;i++)
	{
		for(j=0;j<=m-1;j++)
	        cout<<symmat[i][j]<<" ";
		cout<<endl;
	}

	for(i=0;i<=m-1;i++)
		cout<<evals[i]<<"  ";

    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-1; j >= 0; j--) 
	 {
       cout<<evals[j]<<endl;		//特征值 
	 }

    for(i=0;i<=m-1;i++)
	{
		for(j=0;j<=m-1;j++)
	        cout<<symmat[i][j]<<" ";
		cout<<endl;
	}


     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 = 0; j <= m-1; j++) 
	 {
       for (i = 0; i <= 1; 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 = 0; i <= n-1; i++) 
	 {
      for (j = 0; j <= m-1; j++)
	  {
        interm[j] = data[i][j]; 
	  }   /* data[i][j] will be overwritten */
        for (k = 0; k <= 1; k++) 
		{
          data[i][k] = 0.0;
          for (k2 = 0; k2 <= m-1; k2++)
		  {
            data[i][k] += interm[k2] * symmat[k2][m-k-1]; 
		  }
        }
     }

     printf("\nProjections of row-points on first 3 prin. comps.:\n");
     for (i = 0; i <= n-1; i++)
	 {
       for (j = 0; j <= 1; 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 = 0; j <= m-1; j++) 
	 {
      for (k = 0; k <= m-1; k++) 
	  {
        interm[k] = symmat2[j][k];
	  }  /*symmat2[j][k] will be overwritten*/
        for (i = 0; i <= 1; i++) 
		{
          symmat2[j][i] = 0.0;
          for (k2 = 0; k2 <= m-1; k2++) 
		  {
            symmat2[j][i] += interm[k2] * symmat[k2][m-i-1]; 
		  }
          if (evals[m-i-1] > 0.0005)   /* Guard against zero eigenvalue */
             symmat2[j][i] /= sqrt(evals[m-i-1]);   /* Rescale */
          else
             symmat2[j][i] = 0.0;    /* Standard kludge */
        }
     }

     printf("\nProjections of column-points on first 3 prin. comps.:\n");
     for (j = 0; j <= m-1; j++) 
	 {
       for (k = 0; k <= 1; k++)  
	   {
          printf("%12.4f", symmat2[j][k]);  
	   }
          printf("\n"); 
	 }

    for(i=0;i<=n-1;i++)
		delete[] data[i];
	delete[] data;

    for(i=0;i<=m-1;i++)
		delete[] symmat[i];
	delete[] symmat;

    for(i=0;i<=m-1;i++)
		delete[] symmat2[i];
	delete[] symmat2;

    delete[] evals;

	cin>>t;

}