C/C++实现PCA降维以及故障监测
生活随笔
收集整理的這篇文章主要介紹了
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
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分布
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;//177testdata >> dat.n;//8//新建二維數組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; }部分主函數
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 <<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降維后的維數：" << 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");經過測試，降維結果與Matlabt降維結果相一致，同時進行T^2和Q統計量的閾值與計算值均與Matlab計算結果吻合
程序計算結果與Matlab結果
部分數據對比: 左邊為Matlab計算結果，右側為程序結果，第一幅為T^2統計量對比，第二幅為Q統計量對比 。
通過訓練數據，計算T^2統計量與Q統計量的閾值，通過測試數據計算得到兩個統計量的大小，如上表，經過比較進行故障監測，在此不做詳細介紹。
神經網絡進行故障診斷
BP人工神經網絡
利用PCA降維得到的特征值進行學習。
部分代碼如下：
總結
重要的事說三遍：
沒事別C/C++,老老實實用Python！！！！！！！
沒事別C/C++,老老實實用Python！！！！！！！
沒事別C/C++,老老實實用Python！！！！！！！
總結
以上是生活随笔為你收集整理的C/C++实现PCA降维以及故障监测的全部內容，希望文章能夠幫你解決所遇到的問題。
- 上一篇： 微信小程序 action-sheet组件
- 下一篇： C++ endian转换