📄 pca1a.cpp
字号:
double * Eigen --- 长度为n，为n个特征值
double * EigenVector --- 维数为nxn，返回n阶实对称阵的特征向量组其中，
EigenVector中的j列为与数组Eigen中第j个特征值
对应的特征向量
返回值：
若返回标记小于0，则说明迭代了Iteration次还未求得一个特征值；
若返回标记大于0，则说明程序工作正常，全部特征值由一
维数组Eigen给出，特征向量组由二维数组EigenVector给出
*/
int SymmetricRealMatrix_Eigen( double CovMatrix[], int n,
double Eigen[], double EigenVector[] )
{
int k;
double * subDiagonal;
subDiagonal = ( double * )malloc( sizeof( double ) * n );
Householder_Tri_Symetry_Diagonal( CovMatrix, n, EigenVector, Eigen, subDiagonal );
k = Tri_Symmetry_Diagonal_Eigenvector( n, Eigen, subDiagonal, EigenVector, EPS, Iteration );
free( subDiagonal );
return( k );
}
/*
PCA1: Perform PCA using covariance.
data --- MxN matrix of input data ( M dimensions, N trials )
行数为原始数据维数；每列数据为一个样本
signals --- MxN matrix of projected data
PC --- each column is a PC
V --- Mx1 matrix of variances
row = M dimensions, col = N trials
*/
int pca1( double * data, int row, int col,
double * signals, double * PC, double * V )
{
int x, y, k;
double * mean;
double * data1, * tPC, * tV;
double * covariance, temp;
int rvalue, * no, tp;
if ( row <= 1 || col <= 1 ) return( -1 );
/* subtract off the mean for each dimension */
mean = ( double * )malloc( sizeof( double ) * row );
data1 = ( double * )malloc( sizeof( double ) * row * col );
for ( y = 0; y < row; y++ ) /* calculate mean */
{
mean[y] = 0;
for ( x = 0; x < col; x++ )
mean[y] += data[y*col+x];
}
for ( y = 0; y < row; y++ ) mean[y] = mean[y]/col;
for ( y = 0; y < row; y++ ) /* subtract mean */
for ( x = 0; x < col; x++ )
{
data1[y*col+x] = data[y*col+x] - mean[y];
}
free( mean );
/* calculate the covariance matrix */
covariance = ( double * )malloc( sizeof( double ) * row * row );
for ( y = 0; y < row; y++ )
for ( x = 0; x < row; x++ )
{
temp = 0;
for ( k = 0; k < col; k++ )
temp += data1[y*col+k] * data1[x*col+k];
temp = temp / ( col-1 );
covariance[x*row+y] = temp;
}
/* find the eigenvectors and eigenvalues */
rvalue = SymmetricRealMatrix_Eigen( covariance, row, V, PC );
free( covariance );
/* sort the variances in decreasing order */
no = ( int * )malloc( sizeof( int ) * row );
for ( x = 0; x < row; x++ ) no[x] = x;
for ( x = 0; x < row-1; x++ )
{
temp = V[x];
for ( k = x; k < row; k++ )
if ( temp < V[k] )
{
tp = no[k];
no[k] = no[x];
no[x] = tp;
temp = V[k];
}
}
/* exchange eigen value and vector in decreasing order */
tV = ( double * )malloc( sizeof( double ) * row );
tPC = ( double * )malloc( sizeof( double ) * row * row );
for ( x = 0; x < row; x++ ) tV[x] = V[x];
for ( x = 0; x < row; x++ )
for ( y = 0; y < row; y++ )
tPC[x*row+y] = PC[x*row+y];
for ( x = 0; x < row; x++ )
{
if ( no[x] != x )
{
for ( y = 0; y < row; y++ )
PC[y*row+x] = tPC[y*row+no[x]];
V[x] = tV[no[x]];
}
}
free( no );
free( tV );
free( tPC );
/* project the original data: signals = PC' * data; */
for ( y = 0; y < row; y++ )
for ( x = 0; x < col; x++ )
{
signals[y*col+x] = 0;
for ( k = 0; k < row; k++ )
signals[y*col+x] += PC[k*row+y] * data1[k*col+x];
}
free( data1 );
return( rvalue );
}
/*
投影到主元素向量上
double * newdata --- 要投影到主元素分量上的数据矩阵
int row, col --- 数据矩阵的维数
row（行数）为主元向量维数（即特征数）
col（列数）为采样数
double * PC --- 主元特征向量（列向量）
double * newsignals --- 投影后获得的信号的系数（row*col）
double * ShiftValue --- 如果一个有row个元素的偏移量，如果改值为NULL
则计算均值，并减去均值，否则减去改值
*/
int project2PCA( double * newdata, int row, int col,
double * PC, double * newsignals, double * ShiftValue )
{
int x, y, k;
double * mean, * data1;
/* subtract off the mean for each dimension */
data1 = ( double * )malloc( sizeof( double ) * row * col );
/* if ShiftValue <> NULL */
if ( ShiftValue != NULL )
{
for ( y = 0; y < row; y++ ) /* subtract ShiftValue */
for ( x = 0; x < col; x++ )
{
data1[y*col+x] = newdata[y*col+x] - ShiftValue[y];
}
}
else
{
mean = ( double * )malloc( sizeof( double ) * row );
for ( y = 0; y < row; y++ ) /* calculate mean */
{
mean[y] = 0;
for ( x = 0; x < col; x++ )
mean[y] += newdata[y*col+x];
}
if ( col <= 1 )
for ( y = 0; y < row; y++ ) mean[y] = 0;
else
for ( y = 0; y < row; y++ ) mean[y] = mean[y]/col;
for ( y = 0; y < row; y++ ) /* subtract mean */
for ( x = 0; x < col; x++ )
{
data1[y*col+x] = newdata[y*col+x] - mean[y];
}
free( mean );
}
/* project the original data: signals = PC' * data; */
for ( y = 0; y < row; y++ )
for ( x = 0; x < col; x++ )
{
newsignals[y*col+x] = 0;
for ( k = 0; k < row; k++ )
newsignals[y*col+x] += PC[k*row+y] * data1[k*col+x];
}
free( data1 );
return( 1 );
}
/*
test: main function
*/
int main()
{
double data[37*8] = {
3.0, 3.0, 3.0, 3.0, 3.0, 3.0, 35.0, 45.0,
53.0, 55.0, 58.0, 113.0, 113.0, 86.0, 67.0, 90.0,
3.5, 3.5, 4.0, 4.0, 4.5, 4.5, 46.0, 59.0,
63.0, 58.0, 58.0, 125.0, 126.0, 110.0, 78.0, 97.0,
4.0, 4.0, 4.5, 4.5, 5.0, 5.0, 48.0, 60.0,
68.0, 65.0, 65.0, 123.0, 123.0, 117.0, 87.0, 108.0,
5.0, 5.0, 5.0, 5.5, 5.5, 5.5, 46.0, 63.0,
70.0, 64.0, 63.0, 116.0, 119.0, 115.0, 97.0, 112.0,
6.0, 6.0, 6.0, 6.0, 6.5, 6.5, 51.0, 69.0,
77.0, 70.0, 71.0, 120.0, 122.0, 122.0, 96.0, 123.0,
11.0, 11.0, 11.0, 11.0, 11.0, 11.0, 64.0, 75.0,
81.0, 79.0, 79.0, 112.0, 114.0, 113.0, 98.0, 115.0,
20.0, 20.0, 20.0, 20.0, 20.0, 20.0, 76.0, 86.0,
93.0, 92.0, 91.0, 104.0, 104.5, 107.0, 97.5, 104.0,
30.0, 30.0, 30.0, 30.0, 30.1, 30.2, 84.0, 96.0,
98.0, 99.0, 96.0, 101.0, 102.0, 99.0, 94.0, 99.0,
30.0, 33.4, 36.8, 40.0, 43.0, 45.6, 100.0, 106.0,
106.0, 108.0, 101.0, 99.0, 98.0, 99.0, 95.0, 95.0,
42.0, 44.0, 46.0, 48.0, 50.0, 51.0, 109.0, 111.0,
110.0, 110.0, 103.0, 95.5, 95.5, 95.0, 92.5, 92.0,
60.0, 61.7, 63.5, 65.5, 67.3, 69.2, 122.0, 124.0,
124.0, 121.0, 103.0, 93.2, 92.5, 92.2, 90.0, 90.8,
70.0, 70.1, 70.2, 70.3, 70.4, 70.5, 137.0, 132.0,
134.0, 128.0, 101.0, 91.7, 90.2, 88.8, 87.3, 85.8,
78.0, 77.6, 77.2, 76.8, 76.4, 76.0, 167.0, 159.0,
152.0, 144.0, 103.0, 89.8, 87.7, 85.7, 83.7, 81.8,
98.9, 97.8, 96.7, 95.5, 94.3, 93.2, 183.0, 172.0,
162.0, 152.0, 102.0, 87.5, 85.3, 83.3, 81.3, 79.3,
160.0, 157.0, 155.0, 152.0, 149.0, 147.0, 186.0, 175.0,
165.0, 156.0, 120.0, 87.0, 84.9, 82.8, 80.8, 79.0,
272.0, 266.0, 260.0, 254.0, 248.0, 242.0, 192.0, 182.0,
170.0, 159.0, 131.0, 88.0, 85.8, 83.7, 81.6, 79.6,
382.0, 372.0, 362.0, 352.0, 343.0, 333.0, 205.0, 192.0,
178.0, 166.0, 138.0, 86.2, 84.0, 82.0, 79.8, 77.5,
770.0, 740.0, 710.0, 680.0, 650.0, 618.0, 226.0, 207.0,
195.0, 180.0, 160.0, 82.9, 80.2, 77.7, 75.2, 72.7
};
double newdata[37] = { 3.0, 53.0, 3.5, 63.0, 4.0, 68.0, 5.0, 70.0, 6.0,
77.0, 11.0, 81.0, 20.0, 93.0, 30.0, 98.0, 30.0, 106.0,
42.0, 110.0, 60.0, 124.0, 70.0, 134.0, 78.0, 152.0, 98.9,
162.0, 160.0, 165.0, 272.0, 170.0, 382.0, 178.0, 770.0, 195.0 };
double shiftvalue[37] = { 3.0, 50.0, 3.0, 60.0, 4.0, 60.0, 6.0, 60.0, 4.0,
75.0, 10.0, 80.0, 22.0, 90.0, 20.0, 100.0, 28.0, 100.0,
40.0, 110.0, 50.0, 120.0, 73.0, 130.0, 78.5, 152.0, 90.0,
160.0, 150.0, 160.0, 270.0, 160.0, 372.0, 180.0, 700.0, 185.0 };
int row = 37, col = 8;
double signals[37*8], PC[37*37], V[37], newsignals[37];
int x, y;
pca1( data, row, col, signals, PC, V );
printf( "Project to Principal Component: \n" );
for ( y = 0; y < row; y++ )
{
for ( x = 0; x < col; x++ )
{
printf( "%7.4f ", signals[y*col+x] );
}
printf( "\n" );
}
printf( "Eigen Values (in deceasing order): \n" );
for ( y = 0; y < row; y++ )
printf( "%7.4f ", V[y] );
printf( "\n" );
printf( "Eigen Vectors: \n" );
for ( y = 0; y < row; y++ )
{
for ( x = 0; x < row; x++ )
{
printf( "%7.4f ", PC[y*row+x] );
}
printf( "\n" );
}
printf( "Project a vector to PCA:\n" );
/* project2PCA( newdata, row, 1, PC, newsignals, NULL ); */
project2PCA( newdata, row, 1, PC, newsignals, shiftvalue );
for ( y = 0; y < row; y++ )
printf( "%7.4f ", newsignals[y] );
printf( "\n" );
return( 1 );
}
⌨️ 快捷键说明
复制代码
Ctrl + C
搜索代码
Ctrl + F
全屏模式
F11
切换主题
Ctrl + Shift + D
显示快捷键
?
增大字号
Ctrl + =
减小字号
Ctrl + -