之前使用Matlab進行過關(guān)于PCA的研究，應用于故障診斷與監(jiān)測中，為了方便以后與其他平臺進行耦合，采用C/C++語言實現(xiàn)，參考了很多人編寫的C/C++代碼，也走過一些彎路，針對相關(guān)學習在此記錄下來，有需要的朋友可以參考。

一些降維算法

Principal Component Analysis （PCA）
Linear Discriminant Analysis（LDA）
Locally linear embedding（LLE）
Laplacian Eigenmaps
PCA：PCA算法是一種線性投影技術(shù)，利用降維后使數(shù)據(jù)的方差最大原則保留盡可能多的信息；
KPCA：PCA僅考慮了數(shù)據(jù)的二階統(tǒng)計信息，而沒有利用高階統(tǒng)計信息，忽略了數(shù)據(jù)的非線性相關(guān)性，而KPCA，通過非線性變換將數(shù)據(jù)映射到了高維，在高維空間中進行特征提取，獲得了更好的特征提取性能；
PPCA：PCA沒有將數(shù)據(jù)的概率分布考慮，PPCA對PCA做了概率上的解釋，延伸了PCA算法。
總之：PPCA和KPCA都是針對PCA算法的缺陷，做出了不同方向上的改進。
數(shù)據(jù)降維參考

PCA

PCA（principal components analysis）即主成分分析技術(shù)，又稱主分量分析。主成分分析也稱主分量分析，旨在利用降維的思想，把多指標轉(zhuǎn)化為少數(shù)幾個綜合指標。
在統(tǒng)計學中，主成分分析PCA是一種簡化數(shù)據(jù)集的技術(shù)。它是一個線性變換。這個變換把數(shù)據(jù)變換到一個新的坐標系統(tǒng)中，使得任何數(shù)據(jù)投影的第一大方差在第一個坐標(稱為第一主成分)上，第二大方差在第二個坐標(第二主成分)上，依次類推。主成分分析經(jīng)常用于減少數(shù)據(jù)集的維數(shù)，同時保持數(shù)據(jù)集的對方差貢獻最大的特征。這是通過保留低階主成分，忽略高階主成分做到的。這樣低階成分往往能夠保留住數(shù)據(jù)的最重要方面。但是，這也不是一定的，要視具體應用而定。
PCA是模式識別中常見的特征降維的算法，其大體步驟可以分為以下幾個部分：
（1）原始特征矩陣歸一化處理
（2）求取歸一化處理后特征矩陣的協(xié)方差矩陣
（3）計算協(xié)方差矩陣的特征值及其對應的特征向量
（4）按照特征值從大到小排列特征向量
（5）從大到小，挑選出前K個特征值對應的特征向量組成降維后的特征向量，即為所求。

工具

Boost庫
Boost是為C++語言標準庫提供擴展的一些C++程序庫的總稱。Boost庫是一個可移植、提供源代碼的C++庫，作為標準庫的后備，是C++標準化進程的開發(fā)引擎之一，是為C++語言標準庫提供擴展的一些C++程序庫的總稱。
Boost庫由C++標準委員會庫工作組成員發(fā)起，其中有些內(nèi)容有望成為下一代C++標準庫內(nèi)容。在C++社區(qū)中影響甚大，是不折不扣的“準”標準庫。

Eigen庫
Eigen是一個高層次的C ++庫，有效支持線性代數(shù)，矩陣和矢量運算，數(shù)值分析及其相關(guān)的算法。Eigen是一個開源庫,是一個基于C++模板的線性代數(shù)庫.提供有關(guān)矩陣的的線性代數(shù)運算，解方程等功能。官方的文檔在此，以SLAM十四講代碼閱讀，快速入門。
Eigen庫說明

VS2013
Visual Studio是目前最流行的Windows平臺應用程序的集成開發(fā)環(huán)境。Microsoft Visual Studio（簡稱VS）是美國微軟公司的開發(fā)工具包系列產(chǎn)品。VS是一個基本完整的開發(fā)工具集，它包括了整個軟件生命周期中所需要的大部分工具，如UML工具、代碼管控工具、集成開發(fā)環(huán)境(IDE)等等。
需要在項目屬性中進行Boost和Eigen庫的配置，可以參考下列文檔。
Boost庫配置
Eigen庫配置

部分代碼

由于涉及故障監(jiān)測中的內(nèi)容，需要計算兩個統(tǒng)計量T^2和Q, 其中涉及F分布和正態(tài)分布的計算內(nèi)容，在此不過多介紹。
統(tǒng)計量T^2：設X～Np(μ，∑)，S～Wp(n，∑)，且X與S相互獨立，n≥p，則稱統(tǒng)計量T2=nX’S-1X的分布為非中心Hotelling T2分布，記為T2～T2(p，n，μ)。當μ=0時，稱T2服從(中心)Hotelling T2分布，記為T2(p，n)，由于這一統(tǒng)計量的分布首先由Harold Hotelling提出來的，故稱Hotelling T2分布

#include"pca.h" #include< iostream> #include< fstream> #include<Eigen/Dense> #include<Eigen/Eigenvalues> using namespace Eigen; using namespace std; //*********************************************************** #include<stdlib.h> #include<stdio.h> #include<math.h> 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 <boost\math\distributions\normal.hpp> #include <iostream> #include <random> #include <math.h> //********************************************************** 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 (x<crit) x = crit;if (x>1 - crit) x = 1 - crit;h = 1;while ((h>(crit*fabs(x))) && (h>crit) && (count_max<count_max_limit)){count_max = count_max + 1;if (x>1) 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分布有著廣泛的應用，如在方差分析、回歸方程的顯著性檢驗中都有著重要的地位。
正態(tài)分布使用Boost庫即可，網(wǎng)絡上也可以找到參考，不做介紹。

部分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;//177testdata >> dat.n;//8//新建二維數(shù)組dat.data = new double*[dat.m];for (i = 0; i<dat.m; i++)dat.data[i] = new double[dat.n];for (i = 0; i<dat.m; i++)for (j = 0; j<dat.n; j++)testdata >> 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<columns; i++){s = 0;for (int j = 0; j<rows; j++)s += a[j][i];s = s / rows;ss = 0;for (int j = 0; j<rows; j++)ss += (a[j][i] - s)*(a[j][i] - s);ss = ss / (rows - 1);for (int j = 0; j<rows; j++) a[j][i] = (a[j][i] - s) / sqrt(ss);} } double PCA::product(double *a, double *b) {double sum = 0;for (int i = 0; i<columns; i++)sum += a[i] * b[i];return sum; } double **PCA::matrixproduct(double **a) {int i, j, k;double **c;c = new double*[columns];for (i = 0; i<columns; i++)c[i] = new double[columns];for (i = 0; i<columns; i++)for (j = 0; j<columns; j++){c[i][j] = 0;for (k = 0; k<rows; k++)c[i][j] += a[k][i] * a[k][j];c[i][j] /= (rows - 1);}return c;for (i = 0; i<columns; i++)delete[columns]c[i];delete[columns]c; } void PCA::zhengjiao(double **v) {double **b;double *xx, *yy;int i;xx = new double[columns];yy = new double[columns];b = new double*[columns];for (i = 0; i<columns; i++)b[i] = new double[columns];for (i = 0; i<columns; i++)b[i][0] = v[i][0];for (int j = 1; j<columns; j++)for (i = 0; i<columns; i++){for (int k = 0; k<j; k++){for (int t = 0; t<columns; t++){xx[t] = b[t][k];yy[t] = v[t][j];}b[i][j] = v[i][j] - (product(xx, yy) / product(xx, xx))*b[i][k];}}for (i = 0; i<columns; i++){for (int j = 0; j<columns; j++)xx[j] = b[j][i];yy[i] = sqrt(product(xx, xx));}for (i = 0; i<columns; i++)for (int j = 0; j<columns; j++)v[i][j] = b[i][j] / yy[j];delete[]xx;delete[]yy;for (i = 0; i<columns; i++)delete[columns]b[i];delete[rows]b; }

部分主函數(shù)

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 << "訓練數(shù)據(jù)的行數(shù)為" << m << "，訓練數(shù)據(jù)的列數(shù)為 " << 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 <<VV << endl;//v[i][j] = VV;ofstream foutvv;foutvv.open("vectorresult.txt");foutvv << VV << "\n";foutvv << flush;foutvv.close(); //****************************************************************************FILE*fp2 = fopen("valueresult.txt", "r");if (fp2 == NULL){printf("無法打開文件");}for (i = 0; i<n; i++) //i<500{//j < 21fscanf(fp2, "%lf", &A[i]);}for (i = 0; i<n; i++) //i<500{//j < 21printf("%lf ", A[i]);}fclose(fp2);printf("\n"); //***********************************************cout << "求特征向量矩陣" << endl;system("pause");FILE*fp3 = fopen("vectorresult.txt", "r");if (fp3 == NULL){printf("無法打開文件");}for (i = 0; i<n; i++) //i<500{for (int j = 0; j < n; j++) //j < 21fscanf(fp3, "%lf", &v[i][j]);}fclose(fp3);for (i = 0; i<n; i++) //i<500{for (j = 0; j<n; j++) //j<21{printf("%lf ", v[i][j]);//輸出}printf("\n");} //*****************************************************************************testpca.zhengjiao(v); testpca.selectionsort(A, v); t = testpca.selectcharactor(A, getratio, B); printf("\n");cout << "PCA降維后的維數(shù)：" << t << endl;cout << "排序后提取的特征值" << endl;for (i = 0; i <= t - 1; i++) printf("%13.7e ", A[i]);printf("\n\n");cout << "特征值對應的特征向量" << endl;for (i = 0; i < n; i++) {for (j = 0; j < t; j++)printf("%13.7e ", v[i][j]); printf("\n");} cout << "特征值的累計貢獻率是" << endl;for (i = 0; i < n; i++)cout << B[i] << " ";cout << endl;cout << "當提取效率是" << getratio << "時提取了前" << t << "個分量" << endl; if (t >= 1 && t <= n)Project = testpca.getProject(t, trainx, v); elsecout << "error" << endl; // testpca.saveProject(projectfile, Project, t);cout << "**********************************************" << endl;printf("\n");

經(jīng)過測試，降維結(jié)果與Matlabt降維結(jié)果相一致，同時進行T^2和Q統(tǒng)計量的閾值與計算值均與Matlab計算結(jié)果吻合

程序計算結(jié)果與Matlab結(jié)果

部分數(shù)據(jù)對比: 左邊為Matlab計算結(jié)果，右側(cè)為程序結(jié)果，第一幅為T^2統(tǒng)計量對比，第二幅為Q統(tǒng)計量對比 。


通過訓練數(shù)據(jù)，計算T^2統(tǒng)計量與Q統(tǒng)計量的閾值，通過測試數(shù)據(jù)計算得到兩個統(tǒng)計量的大小，如上表，經(jīng)過比較進行故障監(jiān)測，在此不做詳細介紹。

神經(jīng)網(wǎng)絡進行故障診斷

BP人工神經(jīng)網(wǎng)絡
利用PCA降維得到的特征值進行學習。
部分代碼如下：

#include"stdafx.h" #include <stdio.h> #include <time.h> #include <math.h> #include <stdlib.h> #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; i<Data; i++)for (j = 0; j<In; j++)fscanf(fp1, "%lf", &d_in[i][j]);fclose(fp1);if ((fp2 = fopen("inflag.txt", "r")) == NULL){printf("can not open the inflag.txt\n");exit(0);}for (i = 0; i<Data; i++)for (j = 0; j<Out; j++)fscanf(fp1, "%lf", &d_out[i][j]);fclose(fp2); } //初始化神經(jīng)網(wǎng)絡void initBPNework() {int i, j; for (i = 0; i<In; i++){Minin[i] = Maxin[i] = d_in[0][i];for (j = 0; j<Data; j++){Maxin[i] = Maxin[i]>d_in[j][i] ? Maxin[i] : d_in[j][i];Minin[i] = Minin[i]<d_in[j][i] ? Minin[i] : d_in[j][i];}}//輸出數(shù)據(jù)的最大最小值for (i = 0; i<Out; i++){Minout[i] = Maxout[i] = d_out[0][i];for (j = 0; j<Data; j++){Maxout[i] = Maxout[i]>d_out[j][i] ? Maxout[i] : d_out[j][i];Minout[i] = Minout[i]<d_out[j][i] ? Minout[i] : d_out[j][i];}} for (i = 0; i < In; i++)for (j = 0; j < Data; j++)d_in[j][i] = (d_in[j][i] - Minin[i] + 1) / (Maxin[i] - Minin[i] + 1); for (i = 0; i < Out; i++)for (j = 0; j < Data; j++)d_out[j][i] = (d_out[j][i] - Minout[i] + 1) / (Maxout[i] - Minout[i] + 1); for (i = 0; i < Neuron; ++i)for (j = 0; j < In; ++j){w[i][j] = (rand()*2.0 / RAND_MAX - 1) / 2;dw[i][j] = 0;} for (i = 0; i < Neuron; ++i)for (j = 0; j < Out; ++j){v[j][i] = (rand()*2.0 / RAND_MAX - 1) / 2;dv[j][i] = 0;} }

總結(jié)

重要的事說三遍：

沒事別C/C++,老老實實用Python！！！！！！！

沒事別C/C++,老老實實用Python！！！！！！！

沒事別C/C++,老老實實用Python！！！！！！！

總結(jié)

以上是生活随笔為你收集整理的C/C++实现PCA降维以及故障监测的全部內(nèi)容，希望文章能夠幫你解決所遇到的問題。

如果覺得生活随笔網(wǎng)站內(nèi)容還不錯，歡迎將生活随笔推薦給好友。