C/C++实现PCA降维以及故障监测

之前使用Matlab进行过关于PCA的研究，应用于故障诊断与监测中，为了方便以后与其他平台进行耦合，采用C/C++语言实现，参考了很多人编写的C/C++代码，也走过一些弯路，针对相关学习在此记录下来，有需要的朋友可以参考。

Principal Component Analysis （PCA）
Linear Discriminant Analysis（LDA）
Locally linear embedding（LLE）
Laplacian Eigenmaps
PCA：PCA算法是一种线性投影技术，利用降维后使数据的方差最大原则保留尽可能多的信息；
KPCA：PCA仅考虑了数据的二阶统计信息，而没有利用高阶统计信息，忽略了数据的非线性相关性，而KPCA，通过非线性变换将数据映射到了高维，在高维空间中进行特征提取，获得了更好的特征提取性能；
PPCA：PCA没有将数据的概率分布考虑，PPCA对PCA做了概率上的解释，延伸了PCA算法。
总之：PPCA和KPCA都是针对PCA算法的缺陷，做出了不同方向上的改进。
数据降维参考

PCA（principal components analysis）即主成分分析技术，又称主分量分析。主成分分析也称主分量分析，旨在利用降维的思想，把多指标转化为少数几个综合指标。
在统计学中，主成分分析PCA是一种简化数据集的技术。它是一个线性变换。这个变换把数据变换到一个新的坐标系统中，使得任何数据投影的第一大方差在第一个坐标(称为第一主成分)上，第二大方差在第二个坐标(第二主成分)上，依次类推。主成分分析经常用于减少数据集的维数，同时保持数据集的对方差贡献最大的特征。这是通过保留低阶主成分，忽略高阶主成分做到的。这样低阶成分往往能够保留住数据的最重要方面。但是，这也不是一定的，要视具体应用而定。
PCA是模式识别中常见的特征降维的算法，其大体步骤可以分为以下几个部分：
（1）原始特征矩阵归一化处理
（2）求取归一化处理后特征矩阵的协方差矩阵
（3）计算协方差矩阵的特征值及其对应的特征向量
（4）按照特征值从大到小排列特征向量
（5）从大到小，挑选出前K个特征值对应的特征向量组成降维后的特征向量，即为所求。

Boost库
Boost是为C++语言标准库提供扩展的一些C++程序库的总称。Boost库是一个可移植、提供源代码的C++库，作为标准库的后备，是C++标准化进程的开发引擎之一，是为C++语言标准库提供扩展的一些C++程序库的总称。
Boost库由C++标准委员会库工作组成员发起，其中有些内容有望成为下一代C++标准库内容。在C++社区中影响甚大，是不折不扣的“准”标准库。

Eigen库
Eigen是一个高层次的C ++库，有效支持线性代数，矩阵和矢量运算，数值分析及其相关的算法。Eigen是一个开源库,是一个基于C++模板的线性代数库.提供有关矩阵的的线性代数运算，解方程等功能。官方的文档在此，以SLAM十四讲代码阅读，快速入门。
Eigen库说明

VS2013
Visual Studio是目前最流行的Windows平台应用程序的集成开发环境。Microsoft Visual Studio（简称VS）是美国微软公司的开发工具包系列产品。VS是一个基本完整的开发工具集，它包括了整个软件生命周期中所需要的大部分工具，如UML工具、代码管控工具、集成开发环境(IDE)等等。
需要在项目属性中进行Boost和Eigen库的配置，可以参考下列文档。
Boost库配置
Eigen库配置

由于涉及故障监测中的内容，需要计算两个统计量T^2和Q, 其中涉及F分布和正态分布的计算内容，在此不过多介绍。
统计量T^2：设X～Np(μ，∑)，S～Wp(n，∑)，且X与S相互独立，n≥p，则称统计量T2=nX’S-1X的分布为非中心Hotelling T2分布，记为T2～T2(p，n，μ)。当μ=0时，称T2服从(中心)Hotelling T2分布，记为T2(p，n)，由于这一统计量的分布首先由Harold Hotelling提出来的，故称Hotelling T2分布

#include"pca.h"  
#include< iostream> 
#include< fstream> 
#include
#include
using namespace Eigen;
using namespace std;
//***********************************************************
#include
#include
#include
double FDist(double F, double m, double n);
double betainv(double p, double a, double b);
double betainc(double x, double a, double b);
double beta(double z, double w);
double gamma(double xx);
double beta_cf(double a, double b, double x);
//***************************************************************
#include 
#include 
#include 
#include 
//**********************************************************
double FDist(double F, double m, double n)
{
	double xx, p;

	if (m <= 0 || n <= 0) p = -1;
	else if (F>0)
	{
		xx = F / (F + n / m);
		p = betainc(xx, m / 2, n / 2);
	}
	return(1 - p);
}

double betainv(double p, double a, double b)
{

	int count_max_limit = 100;
	int count_max = 0;
	double x, xnew, y, h, pbeta, logkerna, logkernb;
	double crit = 1.818989403545857e-012; 
	//	float crit=6.4155306e-006;          

	if (p == 0) x = 0;
	if (p == 1) x = 1;

	x = a / (a + b);
	if (x1 - crit) x = 1 - crit;
	h = 1;

	while ((h>(crit*fabs(x))) && (h>crit) && (count_max1) p = 1;
		pbeta = betainc(x, a, b);
		if (pbeta>1) pbeta = 1;

		logkerna = (a - 1)*log(x);
		if ((a == 1) && (x == 0)) logkerna = 0;
		logkernb = (b - 1)*log(1 - x);
		if ((b == 1) && (x == 1)) logkernb = 0;
		y = exp(logkerna + logkernb - log(beta(a, b)));

		h = (pbeta - p) / y;
		xnew = x - h;

		if (xnew <= 0) xnew = x / 10;
		if (xnew >= 1) xnew = 1 - (1 - x) / 10;

		x = xnew;
	}
	return x;
}

F分布参考
F分布C/C++代码
F分布是两个服从卡方分布的独立随机变量各除以其自由度后的比值的抽样分布，是一种非对称分布，且位置不可互换。F分布有着广泛的应用，如在方差分析、回归方程的显著性检验中都有着重要的地位。
正态分布使用Boost库即可，网络上也可以找到参考，不做介绍。

部分pca.h头文件

typedef struct sourcedata   
{
	int m;
	int n;
	double **data;
}SourceData;
class PCA
{
public:
	PCA(int m, int n);   
	SourceData getdata(const  char *file); 
	//double **getdata(const char *file, int &m, int &n);
	void standarddata(double **a); 
	double  product(double *a, double *b); 
	void  swap(double &x, double &y);
	double  **matrixproduct(double **a); 
	void  selectionsort(double *A, double **v);
	void  zhengjiao(double **v); 
	int jcb(double **a, double **v, double eps, int jt);		 
	int selectcharactor(double *A, double getratio, double *B); 
	double  **getProject(int t, double **x, double **v);//计算投影  
	void   saveProject(const char *projectfile, double **project, int t); 
	~PCA(){}
private:
	int  rows;
	int columns;
};
SourceData PCA::getdata(const char *file)
{
	SourceData dat;
	ifstream testdata;
	int i, j;
	testdata.open(file);
	if (!testdata)
	{
		cout << "cannot open" << endl;
	}
	testdata >> dat.m;//177
	testdata >> dat.n;//8
	//新建二维数组
	dat.data = new double*[dat.m];
	for (i = 0; i> dat.data[i][j];
	testdata.close();
	return dat;
}
PCA::PCA(int m, int n)   
{
	columns = n;
	rows = m;
}
void PCA::standarddata(double **a)
{
	double s, ss;
	int i;
	for (i = 0; i 
部分主函数
 
void main()
{
	cout << "-----------------------PCA降维------------------------" << endl;
	int i, j, t;          
	int m, n;    
	int mm, nn;
	double **x, **c, **v, **Project,**trainx;
	double **testx;
	double *A, *B;    
//	double *AA, *BB;
	sourcedata pp;
	sourcedata gg;
	double eps = 0.000001;          
//	double getratio = 0.85;
	double getratio=0; 
	const  char *File = "traindatamn.txt";  
	const  char *Filetest = "testdatamn.txt"; 
//	const char *projectfile = "pcapdata.txt";    
	PCA pca(2, 3);    
	pp = pca.getdata(File);   
	gg = pca.getdata(Filetest);
	x = pp.data;
	trainx = pp.data;
	m = pp.m;
	n = pp.n;
	testx = gg.data;
	mm = gg.m;
	nn = gg.n;
	cout << "训练数据的行数为" << m << "，训练数据的列数为 " << n << endl;
	A = new double[n];                       
	B = new double[n];                      
	v = new double*[n];
	for (i = 0; i < n; i++)
	v[i] = new double[n];
	PCA  testpca(m, n);          
//*********************************************************************
	printf("请输入特征值提取率：");
	scanf("%lf", &getratio); 
	printf("n");
		MatrixXd VV = es.pseudoEigenvectors();
//	cout << "求对角线特征值矩阵" << endl;
//	cout << D << endl;
//	printf("n");
//	cout << "求特征向量矩阵" << endl;
//	cout <= 1 && t <= n)
		Project = testpca.getProject(t, trainx, v); 
	else
		cout << "error" << endl;
//	testpca.saveProject(projectfile, Project, t);
   cout << "**********************************************" << endl;
   printf("n");
 
经过测试，降维结果与Matlabt降维结果相一致，同时进行T^2和Q统计量的阈值与计算值均与Matlab计算结果吻合
 
程序计算结果与Matlab结果 
部分数据对比: 左边为Matlab计算结果，右侧为程序结果，第一幅为T^2统计量对比，第二幅为Q统计量对比 。
 
 
 通过训练数据，计算T^2统计量与Q统计量的阈值，通过测试数据计算得到两个统计量的大小，如上表，经过比较进行故障监测，在此不做详细介绍。
 
神经网络进行故障诊断 
BP人工神经网络
 利用PCA降维得到的特征值进行学习。
 部分代码如下：
 
#include"stdafx.h"
#include  
#include  
#include  
#include  
#define Data  100      
#define In 4           
#define Out 1         
#define Neuron 8       
#define TrainC 10000    
#define A  0.2        
#define B  0.4 
#define a  0.2
#define b  0.3
double d_in[Data][In], d_out[Data][Out];                
double w[Neuron][In], o[Neuron], v[Out][Neuron];        
double Maxin[In], Minin[In], Maxout[Out], Minout[Out];
double OutputData[Out];                               
double dv[Out][Neuron], dw[Neuron][In];                
double e;

int  main(int argc, char const *argv[])
{
	                      
	readData();                         
	initBPNework();                    
	trainNetwork();                   
	printf("   11 0.68 0.38 0.18    "  " %lf n", result(11,0.68,0.38,0.18));   
	printf("   10 1.729 0.4 0.3   "  " %lf n", result(10,1.729, 0.4, 0.3)); 
	printf("   6.4 3.59 1.19 0.67    "  " %lf n", result(6.4,3.59,1.19,0.67)); 
	printf("   6.6 3.44 1.23 0.657    "  " %lf n", result(6.6,3.44,1.23,0.657)); 
	printf("   11.0878130000 1.0340320000 0.4226859600 0.2784231900 "  " %lf n", result(11.0878130000,1.0340320000,0.4226859600,0.2784231900)); 
	  
	system("pause");

	return 0;

	//  system("pause");
}
void readData()
{
	FILE *fp1, *fp2;
	int i, j;
	if ((fp1 = fopen("indata.txt", "r")) == NULL)
	{
		printf("can not open the indata.txt n");
		exit(0);
	}
	for (i = 0; id_in[j][i] ? Maxin[i] : d_in[j][i];
			Minin[i] = Minin[i]d_out[j][i] ? Maxout[i] : d_out[j][i];
			Minout[i] = Minout[i] 
总结 
重要的事说三遍：
 
没事别C/C++,老老实实用Python！！！！！！！ 
没事别C/C++,老老实实用Python！！！！！！！ 
没事别C/C++,老老实实用Python！！！！！！！