第9讲--多元线性回归--主成分回归

1、多元线性回归主成分及其回归第9讲多元线性回归解决的问题系数矩阵系数矩阵 Y=XA建模：求解建模：求解回归系数回归系数A，该过程称为建，该过程称为建模模预报：在预报：在A已知时，对于新测已知时，对于新测Xnew，预报预报Ynew，称为称为预报预报例子某保健品含片产品，说明书标明：由营养某保健品含片产品，说明书标明：由营养物质物质A、B、C组成，产品标注中写出了每组成，产品标注中写出了每片中片中A、B、C物质的含量。问，如何认定物质的含量。问，如何认定？配置配置A、B、C的一组溶液，建立浓度与光的一组溶液，建立浓度与光吸收的关系。既建模求回归系数吸收的关系。既建模求回归系数将药片配置成溶液，测吸光

2、，利用上面的将药片配置成溶液，测吸光，利用上面的模型，预报浓度。模型，预报浓度。建模公式推导Y=XAXtY=XtXA(XtX)-1XtY=AE:学校教学pythonX.txtE:学校教学pythonY.txt问题求解的关键步骤是什么？问题求解的关键步骤是什么？方程数与未知数的关系设有规律上符合如下方程的一设有规律上符合如下方程的一 组实验数据组实验数据 y= ax+b通过实验，不断变更通过实验，不断变更x，测得对应的，测得对应的y求求a,b的值，需要几组这样的数据？的值，需要几组这样的数据？唯一解唯一解最小二乘解最小二乘解y1y2ynx1 1x2 1 1xn 1ab=矩阵形式矩阵形式XtX是是

3、2*2的矩阵的矩阵方程数与未知数的关系设有规律上符合如下方程的一设有规律上符合如下方程的一 组实验数据组实验数据 y= 1.2 x1 + 0.9 x2 + 3.3 x3通过实验，不断变更通过实验，不断变更x1、x2、x3，测得对应的，测得对应的y需要几组这样的数据？需要几组这样的数据？唯一解唯一解最小二乘解最小二乘解方程数小于未知数，一定无解吗y= 1.2 x1 + 0.9 x2 + 3.3 x3当当X1,X2,X3存在线性相关时，问题会怎样？存在线性相关时，问题会怎样？如果如果x个数很多，样本打不到要求，怎么办个数很多，样本打不到要求，怎么办？现实中存在这样的问题吗不同浓度成分相同的溶液，在

4、不同波长不同浓度成分相同的溶液，在不同波长x1、x2下下的吸光值的比值，溶液浓度变化，比值不变。的吸光值的比值，溶液浓度变化，比值不变。既既X1和和X2之间是线性相关的。之间是线性相关的。怎样知道变量之间有相关性？答案：通过线性变化答案：通过线性变化主成分算法能解决这类问题主成分算法能解决这类问题死计算：检查死计算：检查XtX有没有逆，没逆，则线性相关有没有逆，没逆，则线性相关10主成份分析主成份分析 PCA Principle Component Analysis能有效的提取测量数据的有用信息能有效的提取测量数据的有用信息解决变量之间的相关性问题解决变量之间的相关性问题有效去除误差，建立有效

5、的模型有效去除误差，建立有效的模型11PCA分解算法原理分解算法原理采用非线性迭代偏最小二乘法采用非线性迭代偏最小二乘法(Nonlinear Iterative Partial Least Squares, NIPALS)方法方法分解量测矩阵分解量测矩阵S S = T Pt + E =tipi + E主成分示例主成分示例12方差最大方向方差最大方向NIPALS算法每次只求一个主成分，目前最大散差方向算法每次只求一个主成分，目前最大散差方向仪器的信噪比仪器的信噪比仪器测量时，信号强度要远远大于噪声仪器测量时，信号强度要远远大于噪声信号的数据的方差要远远大于噪声的方差信号的数据的方差要远远大于噪声

6、的方差所以，所以，PCA可以区别噪声可以区别噪声样例样例x0.91.10.80.87 22.21.92.1y1.21.00.92 1.11.81 1.91.72.5原数据图PCA后15通过特征值比值判断有效变量数通过特征值比值判断有效变量数在在i/ i+i，应该达到最大应该达到最大值值根据根据i值，取值，取T和和P的前的前i列，即可扔掉噪声列，即可扔掉噪声16主成分回归主成分回归PCRPrinciple Component Regression是多元线性回归！原来 Y=XA 现在 Y=TAT为X的主成分得分，即X经PCA分解后的得分因为T只是X的线性组合，提取了线性相关的部分，且只取前i列，所

7、以模型稳定，去掉噪声numpy中主成分中主成分分解分解SVD分解分解实矩阵的实矩阵的SVD(Singular Value Decomposition，奇，奇异值分解异值分解 )分解分解：分解结果：分解结果：A=USV其中其中S是对角矩阵是对角矩阵numpy中主成分中主成分分解分解-SVD程序代码：程序代码： B = np.linalg.svd(A,full_matrices=False)full_matrices=False一定要写，否则会按复数一定要写，否则会按复数分解分解分解结果：U=B0lamda=B1V = B2Lamda是所有的特征值，可以计算相邻比值，决定主成分，它不是一个矩阵实例

8、光谱矩阵的SVD分解数据：E:学校教学教改项目教材数据S-093790.txt是一个16*6的矩阵看看能求解个特征值？16个？ 6个？96个？实例光谱矩阵的SVD分解data=np.mafromtxt(E:学校教学教改项目教材数据S-093790.txt)data=data.dataB = np.linalg.svd(data,full_matrices=False) B1array( 5.48250094e+00, 1.10440342e+00, 3.27012276e-01, 3.23153080e-03, 2.19720845e-03, 1.11546885e-03)实例光谱矩阵的SVD

9、分解 ld = B1 for i in range(len(ld)-1):temp = ldi/ldi+1print (temp)4.964219451633.37725370576101.1942313561.470743841531.96976226926矩阵中有3个有效特征值根据有效特征值，设定PCA的得分和载荷实例光谱矩阵的SVD分解根据主成分，规划得分U和载荷矩阵PSVD ：X=USVPCA：X=TPtT=US ，P=Vt i=len(lamda) S=np.zeros (i,i) S:i,:i=np.diag (lamda) T = np.dot (U,S)V=V.TP = VT

10、= T:,:kP = P:,:k可否编写PCA类传递矩阵给类求得T、P矩阵，特征值比值列表根据特征值比值，规划T和PPCA类import numpy as npclass PCA: def _init_(self, A): self.A=A def SVDdecompose(self): B = np.linalg.svd(self.A,full_matrices=False) U=B0 lamda=B1 V = B2 i=len(lamda) S=np.zeros (i,i) S:i,:i=np.diag (lamda)PCA类 self.T = np.dot (U,S) V=V.T sel

11、f.P = V compare= for i in range(len(lamda)-1): temp = lamdai/lamdai+1 compare.append(temp) return U,S,V,compare def PCAdecompose(self,k): T = self.T:,:k P = self.P:,:k return T,PPCA类调用p应该先调用decompose方法，根据返回的特征之比值，确定主成分p再调用PCAdecompose方法，设定得分和载荷矩阵主成分回归原来 Y=XA现在 Y=TATtY=TtTATtT-1TtY=A求得最终的回归系数：主成分 X=T

12、Pt因为P是正交矩阵，所以 T=XPY=TA=XPA=XAnew主成分回归数据E:学校教学pythonS-093843.txtE:学校教学pythonC-093843.txt求解方程 C=SAS是6*16的矩阵 所以StS的逆不存在S是光谱矩阵，光谱的不同波长间线性相关，所以可以用PCR程序代码建模过程S=np.mafromtxt(“E:学校教学pythonS-093843.txt )S=S.dataC=np.mafromtxt(“E:学校教学pythonC-093843.txt )C=C.dataB = np.linalg.svd(data,full_matrices=False)U=B0l

13、amda=B1i=len(lamda)S=np.zeros (i,i)S:i,:i=np.diag (lamda)T = np.dot (U,S)V=B2程序代码建模过程P=V.Tfor i in range(len(lmada)-1):temp = lamdai/lamdai+1print (temp)k=int(input(“主成分数为:”)T = T:,:kP = P:,:kTtT=np.dot(T.T,T)inv = np.linalg.inv(TtT)A=np.dot(inv, T.T)Alast=np.dot(P,A)程序代码预报对新测定Snew:C=np.dot(Snew, Al

14、ast)扩展-能否用MLR、PCA类PCR类传递X，Y给PCRPCR内，以X调用PCA，确定主成分数根据确定的主成分数，确定T、P，以T，Y建模，并结合P确定回归系数建立预报方法。扩展-能否用MLR、PCA类PCR类import numpy as npfrom PCA import PCAfrom MLR import MLRclass PCR: def _init_(self,X,Y): self.X=X self.Y=Y def confirmPCs(self): pca=PCA(self.X) U,S,V,compare=pca.SVDdecompose() return compare扩展-能否用MLR、PCA类def model(self,PCs): pca=PCA(self.X) U,S,V,compare=pca.SVDdecompose() T,P=pca.PCAdecompose(PCs) mlr=MLR(T,self.Y) mlr.modelling() self.A=np.dot(P,mlr.A) def predict(self,Xnew): ans = np.dot(Xnew,self.A) return ans调用函数S=np.mafromt