PCA故障诊断的Python实现

# # *****************加载程序所需的函数工具箱 ********************* import pandas as pd import math import numpy import numbers import numpy as np from numpy import * import numpy as transpose from numpy.linalg import * from pylab import * from scipy.stats import f import matplotlib.pyplot as plt import os from pandas import DataFrame # *****************使用pandas方法读取样本数据功能模块（开始）********************* train_data = pd.read_csv('train.csv',header=None) test_data = pd.read_csv('test.csv',header=None) # *****************使用pandas方法读取样本数据功能模块（结束）********************* m = train_data.shape[1];#获取数据表格的列数 n = train_data.shape[0];#获取数据表格的行数 # ******************数据标准化处理（开始）********************* S_mean = np.mean(train_data, axis=0)#健康数据矩阵的列均值 S_mean = np.array(S_mean)#健康数据的列均值，narry数据类型 S_var = np.std(train_data, ddof=1);#健康数据矩阵的列方差,默认ddof=0表示对正态分布变量的方差的最大似然估计，ddof=1提供了对无限总体样本的方差的无偏估计（与Matlab一致） S_var[S_var == 0.0] = 0.0000000000000001#将集合S_var中的0替换为0.0000000000000001 S_var = np.array(S_var)#健康数据的列方差，narry数据类型 train_data -= S_mean #求取矩阵X的均值 train_data /= S_var #求取矩阵X的方差 train_data = np.where(train_data < 4.0e+11,train_data,0.0) #把标准化后的矩阵X中的0替换为0.0000000000000001 X_new = train_data;#求得标准化处理后的矩阵X_new # ******************求矩阵Y的协方差矩阵Z ********************* X_new = np.transpose(X_new);#对矩阵进行转秩操作 Z = np.dot(X_new,train_data/(n-1))#求取协方差矩阵Z # ******************计算协方差矩阵Z的特征值和特征向量 ********************* a,b = numpy.linalg.eig(Z) ##特征值赋值给a，对应特征向量赋值给b lambda1 = sorted(a, reverse = True) #特征值从大到小排序 # ******************主元个数选取 ********************* totalvar = 0; #累计贡献率，初值0 for i in range(m): totalvar = totalvar + lambda1[i]/sum(a)#累计贡献率，初值0 if totalvar >= 0.85: k = i+1 #确定主元个数 break#跳出for循环 PCnum = k; #选取的主元个数 PC = np.eye(m,k);#定义一个矩阵，用于存放选取主元的特征向量 for j in range(k): wt = a.tolist().index(lambda1[j])#查找排序完成的第j个特征值在没排序特征值里的位置。 PC[:,j:j+1] = b[:,wt:wt+1];#提取的特征值对应的特征向量 print(PC) # ******************根据建模数据求取 T2 阈值限 ********************* # ******************置信度 = (1-a)% =（1-0.05）%=95% ************* F =f.ppf(1 - 0.05, k, n-1) #F分布临界值 T2 = k*(n-1)*F/(n-k) #T2求取 # ****************** 健康数据的 SPE 阈值限求解 ********************* ST1 = 0 #对应SPE公式中的角1初值 ST2 = 0 #对应SPE公式中的角2初值 ST3 = 0 #对应SPE公式中的角3初值 for i in range(k-1,m): ST1 = ST1 + lambda1[i] #对应SPE公式中的角1 ST2 = ST2 + lambda1[i]*lambda1[i]#对应SPE公式中的角2 ST3 = ST3 + lambda1[i]*lambda1[i]*lambda1[i]#对应SPE公式中的角3 h0 = 1-2*ST1*ST3/(3*pow(ST2,2)) Ca = 1.6449 SPE = ST1*pow(Ca*pow(2*ST2*pow(h0,2),0.5)/ST1+1+ST2*h0*(h0-1)/pow(ST1,2),1/h0)#健康数据SPE计算 # ******************测试样本数据********************* m1 = test_data.shape[1];#获取数据表格的列数 n1 = test_data.shape[0];#获取数据表格的行数 test_data = np.array(test_data) # 将DataFrame数据烈性转化为ndarray类型，使得数据矩阵与Matlab操作一样。 I = np.eye(m)#产生m*m的单位矩阵 PC1 = np.transpose(PC)#PC的转秩 SPEa = np.arange(n1).reshape(1, n1)#定义测试数据的SPE矩阵,为正数矩阵 SPEa = np.double(SPEa)#将正数矩阵，转化为双精度数据矩阵 TT2a = np.arange(n1).reshape(1, n1)#定义测试数据的T2矩阵,为正数矩阵 TT2a = np.double(TT2a)#将正数矩阵，转化为双精度数据矩阵 DL = np.diag(lambda1[0:k])#特征值组成的对角矩阵 DLi = np.linalg.inv(DL)#特征值组成的对角矩阵的逆矩阵 # ******************绘制结果 ********************* mpl.rcParams['font.sans-serif'] = ['SimHei']#在图形中显示汉字 for i in range(n1): xnew = (test_data[i,:] - S_mean)/S_var;#对应 Matlab程序：xnew=(Data2(i,1:m)-S_mean)./S_var; #以下是实现Matlb程序： err(1,i)=xnew*(eye(14)-PC*PC')*xnew'; xnew1 = np.transpose(xnew)#xnew的转秩 PC1 = np.transpose(PC)#PC的转秩 XPC = np.dot(xnew,PC)#矩阵xnew与PC相乘 XPCPC1 = np.dot(XPC,PC1)#矩阵XPC与PC1相乘 XXPCPC1 = xnew - XPCPC1#矩阵xnew减去XPCPC1 SPEa[0,i] = np.dot(XXPCPC1,XXPCPC1)#矩阵XXPCPC1与XXPCPC1相乘 XPi = np.dot(XPC,DLi)#矩阵XPC与DLi相乘 XPiP = np.dot(XPi,PC1)#矩阵XPi与PC1相乘 TT2a[0,i] = np.dot(XPiP,xnew1)#矩阵XPiP与xnew1相乘 Sampling = r_[0.:n1]#产生的序列值式0到n1 SPE1 = SPE*ones((1,n1))#产生SPE数值相同的矩阵 T21 = T2*ones((1,n1))#产生T2数值相同的矩阵 ########################################################################################### figure(1) #画的第一张图 plot(Sampling, TT2a[0,:], '*-', Sampling,T21[0,:], 'r-')#绘制出测试数据SPEa的数据集合，和健康数据训练得到的SPE阈值限 xlabel('样本个数')#给X轴加标注 ylabel('T^2')#给Y轴加标注 legend(['T^2 检测值','T^2 阈值限'])#为绘制出的图形线条添加标签注明 title("T^2统计量")#绘制的图形主题为“SPE统计量” # show()#显示绘制的图形 figure(2) plot(Sampling, SPEa[0,:],'*-', Sampling,SPE1[0,:], 'r-')#绘制出测试数据TT2a的数据集合，和健康数据训练得到的T2阈值限 xlabel('样本个数')#给X轴加标注 ylabel('SPE')#给Y轴加标注 legend(['SPE 检测值','SPE 阈值限'])#为绘制出的图形线条添加标签注明 title("SPE统计量")#绘制的图形主题为“SPE统计量” show()#显示绘制的图形 #########构建循环输出###################################################################### #循环对象TT2a,SPEa,循环基线T2,SPE sum=0 for i in range(n1):#对测试样本个数进行循环 if ((TT2a[0,i]<=T2)&(SPEa[0,i]<=SPE)):#判断各个值是否小于阈值线 TT2a[0,i]=0#将小于阈值线的样本点位置上的数置为0 SPEa[0,i]=0#将小于阈值线的样本点位置上的数置为0 else: TT2a[0, i] = 1#将小于阈值线的样本点位置上的数置为1 SPEa[0, i] = 1#将小于阈值线的样本点位置上的数置为1 sum+=1 # print(i)#输出有故障的样本点 print(sum) ########################################################################################### d1=DataFrame(TT2a.T) d1['label']=d1[0] d1.drop(0,axis=1,inplace=True) d1.to_csv('label.csv',index=False) print(d1.sum()) print(SPEa) ##########################################################################################