main.cpp

PRINCIPAL COMPONENT ANALYSIS,using C++ language

#include "MyMatrix.h"
#include <vector>

//principal components analysis
int main()
{
	int sampleNum=9,variableNum=6;
	int i,j;
	const double contributionRate=0.85;
	double originalSamples[MAXSIZE][MAXSIZE]={
		{11.5,9.0,7.5,9.1,11.6,13.0,11.6,10.7,11.1},
		{95.3,97.7,110.7,89.0,88.0,87.7,79.7,119.3,87.7},
		{26.4,30.8,39.7,35.4,29.3,24.6,25.6,29.9,32.2},
		{39.2,46.8,39.1,35.3,37.0,44.8,43.7,38.8,35.6},
		{109,107,108,110,109,109,109,108},
		{6.0/2.0,6.0/6.0,6.0/3.0,6.0/2.0,6.0/7.0,6.0/5.0,6.0/5.0,6.0/3.0}
	};
	double a[MAXSIZE][MAXSIZE]={
			{1.0,-0.4933,-0.8982,0.1448,0.8127,0.4397},
			{-0.4933,1.0,0.4369,-0.0853,-0.2979,-0.1263},
			{-0.8982,0.4369,1.0,-0.4709,-0.6883,-0.4925},
			{0.1448,-0.0853,-0.4709,1.0,-0.1653,0.5684},
			{0.8127,-0.2979,-0.6883,-0.1653,1.0,0.3823},
			{0.4397,-0.1263,-0.4925,0.5684,0.3823,1.0}
	};
	double s[MAXSIZE][MAXSIZE]={
		{2,-1,0},
		{-1,2,-1},
		{0,-1,2}
	};
	MyMatrix Sa(s,3,3),Z;
	Sa.JacobiMethod(Z).Display();
	Z.Display();

	MyMatrix originalSampleMatrix(originalSamples,variableNum,sampleNum);

	//计算相关阵时无需对样本进行标准化
	MyMatrix corelationMatrix=originalSampleMatrix.GetCorelationMatrix();
//	MyMatrix corelationMatrix(a,variableNum,variableNum);
	cout<<"相关阵："<<endl;
	corelationMatrix.Display();

	//计算主成分得分值时需要使用标准化后的值
	MyMatrix standarlizedSampleMatrix=originalSampleMatrix.StandarlizeVectors();

	//使用雅克比法计算特征值及对应特征向量
	MyMatrix eigenValueMatrix,eigenVectors;
	eigenValueMatrix=corelationMatrix.JacobiMethod(eigenVectors);

	//提取特征值至数组eigenValues
	double eigenValues[MAXSIZE];
	for(i=0;i<variableNum;i++){
		eigenValues[i]=eigenValueMatrix.GetElement(i,i);
	}

	//从大到小对特征值进行冒泡排序，同时相应调整特征向量顺序
	for(i=0;i<eigenValueMatrix.GetCol();i++){
		for(j=0;j<eigenValueMatrix.GetCol();j++){
			if(eigenValues[i]>eigenValues[j]){
				double t=eigenValues[i];
				eigenValues[i]=eigenValues[j];
				eigenValues[j]=t;
				eigenVectors.SwapColumn(i,j);
			}
		}
	}

	cout<<endl<<"样本相关阵的特征值及特征向量："<<endl;
	MyMatrix t(NULL,1,variableNum);
	for(i=0;i<variableNum;i++)
		t.SetElement(0,i,eigenValues[i]);
	t.Display();
	cout<<endl;
	eigenVectors.Display();

	double sum=0,rate=0;
	for(i=0;i<variableNum;i++){
		sum+=eigenValues[i];
	}

	for(i=0;i<variableNum;i++){
		rate+=eigenValues[i]/sum;
		if(rate>contributionRate){
			cout<<"取前"<<(i+1)<<"个主成分即可"<<endl;
			break;
		}
	}

	//计算得分值，注意样本的指标值是用标准化后的值
	MyMatrix scoreMatrix=eigenVectors*standarlizedSampleMatrix;
	scoreMatrix=scoreMatrix.TransposeMatrix();
	cout<<"各样本的主成分得分值："<<endl;
	scoreMatrix.Display();

	//MyMatrix A(a,6,6);
	//cout<<"Using Jacobi method to get all eigenvalues of matrix C ,"<<endl
	//	<<"diagonal elements of the showed matrix are eigenvalues"<<endl;
	//MyMatrix eigenVectors;
	//MyMatrix C=A.JacobiMethod(eigenVectors);
	//C.Display();
	////eigenVectors.Display();

	//double eigenValues[MAXSIZE];
	//for(int i=0;i<C.GetCol();i++){
	//	eigenValues[i]=C.GetElement(i,i);
	//}

	//for(int i=0;i<C.GetCol();i++){
	//	for(int j=0;j<C.GetCol();j++){
	//		if(eigenValues[i]>eigenValues[j]){
	//			double t=eigenValues[i];
	//			eigenValues[i]=eigenValues[j];
	//			eigenValues[j]=t;
	//			eigenVectors.SwapColumn(i,j);
	//		}
	//	}
	//}

	//MyMatrix t(NULL,1,C.GetCol());
	//for(int i=0;i<C.GetCol();i++)
	//	t.SetElement(0,i,eigenValues[i]);

	//t.Display();
	//eigenVectors.Display();
	return 0;
}