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

多元线性回归主成分及其回归第9讲多元线性回归解决的问题系数矩阵Y=XA建模：求解回归系数A，该过程称为建模预报：在A已知时，对于新测Xnew，预报Ynew，称为预报例子某保健品含片产品，说明书标明：由营养物质A、B、C组成，产品标注中写出了每片中A、B、C物质的含量。问，如何认定？配置A、B、C的一组溶液，建立浓度与光吸收的关系。既建模求回归系数将药片配置成溶液，测吸光，利用上面的模型，预报浓度。建模公式推导Y=XAXtY=XtXA(XtX)-1XtY=AE:\学校教学\python\X.txtE:\学校教学\python\Y.txt问题求解的关键步骤是什么？方程数与未知数的关系设有规律上符合如下方程的一组实验数据y=ax+b通过实验，不断变更x，测得对应的y求a,b的值，需要几组这样的数据？唯一解最小二乘解y1y2…ynx11x21…1xn1ab=矩阵形式XtX是2*2的矩阵方程数与未知数的关系设有规律上符合如下方程的一组实验数据y=1.2x1+0.9x2+3.3x3通过实验，不断变更x1、x2、x3，测得对应的y需要几组这样的数据？唯一解最小二乘解方程数小于未知数，一定无解吗y=1.2x1+0.9x2+3.3x3当X1,X2,X3存在线性相关时，问题会怎样？如果x个数很多，样本打不到要求，怎么办？现实中存在这样的问题吗不同浓度成分相同的溶液，在不同波长x1、x2下的吸光值的比值，溶液浓度变化，比值不变。既X1和X2之间是线性相关的。怎样知道变量之间有相关性？答案：通过线性变化主成分算法能解决这类问题死计算：检查XtX有没有逆，没逆，则线性相关10主成份分析PCAPrincipleComponentAnalysis能有效的提取测量数据的有用信息解决变量之间的相关性问题有效去除误差，建立有效的模型11PCA分解算法原理采用非线性迭代偏最小二乘法(NonlinearIterativePartialLeastSquares,NIPALS)方法分解量测矩阵SS=TPt+E=Σtipi+ET得分矩阵特征值方程Ax=λxP载荷矩阵T和P都是列正交矩阵T的第i列ti的模，就是第i个特征值λiE为残差矩阵，对应噪声每个主成分就是T和P的对应列主成分示例12方差最大方向NIPALS算法每次只求一个主成分，目前最大散差方向仪器的信噪比仪器测量时，信号强度要远远大于噪声信号的数据的方差要远远大于噪声的方差所以，PCA可以区别噪声样例x0.91.10.80.8722.21.92.1y1.21.00.921.11.811.91.72.5t11.4861.4851.2161.3932.6942.8982.5453.253t2-0.2080.075-0.081-0.1580.1420.2210.149-0.273原数据图PCA后15通过特征值比值判断有效变量数在λi/λi+i，应该达到最大值根据i值，取T和P的前i列，即可扔掉噪声16主成分回归PCRPrincipleComponentRegression是多元线性回归！原来Y=XA现在Y=TAT为X的主成分得分，即X经PCA分解后的得分因为T只是X的线性组合，提取了线性相关的部分，且只取前i列，所以模型稳定，去掉噪声numpy中主成分分解—SVD分解实矩阵的SVD(SingularValueDecomposition，奇异值分解)分解：分解结果：A=USV其中S是对角矩阵numpy中主成分分解---SVD程序代码：B=np.linalg.svd(A,full_matrices=False)full_matrices=False一定要写，否则会按复数分解分解结果：U=B[0]lamda=B[1]V=B[2]Lamda是所有的特征值，可以计算相邻比值，决定主成分，它不是一个矩阵实例—光谱矩阵的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)B[1]array([5.48250094e+00,1.10440342e+00,3.27012276e-01,3.23153080e-03,2.19720845e-03,1.11546885e-03])实例—光谱矩阵的SVD分解ld=B[1]foriinrange(len(ld)-1):temp=ld[i]/ld[i+1]print(temp)4.964219451633.37725370576101.1942313561.470743841531.96976226926矩阵中有3个有效特征值根据有效特征值，设定PCA的得分和载荷实例—光谱矩阵的SVD分解根据主成分，规划得分U和载荷矩阵PSVD：X=USVPCA：X=TPtT=US，P=Vti=len(lamda)S=np.zeros((i,i))S[:i,:i]=np.diag(lamda)T=np.dot(U,S)V=V.TP=VT=T[:,:k]P=P[:,:k]可否编写PCA类传递矩阵给类求得T、P矩阵，特征值比值列表根据特征值比值，规划T和PPCA类importnumpyasnpclassPCA:def__init__(self,A):self.A=AdefSVDdecompose(self):B=np.linalg.svd(self.A,full_matrices=False)U=B[0]lamda=B[1]V=B[2]i=len(lamda)S=np.zeros((i,i))S[:i,:i]=np.diag(lamda)PCA类self.T=np.dot(U,S)V=V.Tself.P=Vcompare=[]foriinrange(len(lamda)-1):temp=lamda[i]/lamda[i+1]compare.append(temp)returnU,S,V,comparedefPCAdecompose(self,k):T=self.T[:,:k]P=self.P[:,:k]returnT,PPCA类调用应该先调用decompose方法，根据返回的特征之比值，确定主成分再调用PCAdecompose方法，设定得分和载荷矩阵主成分回归原来Y=XA现在Y=TATtY=TtTATtT-1TtY=A求得最终的回归系数：主成分X=TPt因为P是正交矩阵，所以T=XPY=TA=XPA=XAnew主成分回归数据E:\学校教学\python\S-093843.txtE:\学校教学\python\C-093843.txt求解方程C=SAS是6*16的矩阵所以StS的逆不存在S是光谱矩阵，光谱的不同波长间线性相关，所以可以用PCR程序代码—建模过程S=np.mafromtxt(“E:\\学校教学\\python\\S-093843.txt)S=S.dataC=np.mafromtxt(“E:\\学校教学\\python\\C-093843.txt)C=C.dataB=np.linalg.svd(data,full_matrices=False)U=B[0]lamda=B[1]i=len(lamda)S=np.zeros((i,i))S[:i,:i]=np.diag(lamda)T=np.dot(U,S)V=B[2]程序代码—建模过程P=V.Tforiinrange(len(lmada)-1):temp=lamda[i]/lamda[i+1]print(temp)k=int(input(“主成分数为:”))T=T[:,:k]P=P[:,:k]TtT=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,Alast)扩展--能否用MLR、PCA类PCR类•传递X，Y给PCR•PCR内，以X调用PCA，确定主成分数•根据确定的主成分数，确定T、P，以T，Y建模，并结合P确定回归系数•建立预报方法。扩展--能否用MLR、PCA类PCR类importnumpyasnpfromPCAimportPCAfromMLRimportMLRclassPCR:def__init__(self,X,Y):self.X=Xself.Y=YdefconfirmPCs(self):pca=PCA(self.X)U,S,V,compare=pca.SVDdecompose()returncompare扩展----能否用MLR、PCA类defmodel(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)defpredict(self,Xnew):ans=np.dot(Xnew,self.A)returnans调用函数S=np.mafromtxt(E:\\学校教学\\python\\S-093843.txt)S=S.dataS=S.T要一行一个样本，所以转置C=np.mafromtxt(E:\\学校教学\\python\\C-093843.txt)C=C.dataC=C.Tpcr=PCR(S,C)print(相邻特征值比值)compare=pcr.confirmPCs()print(compare)k=int(input(确定主成分数：))pcr.model(k)ans=pcr.predict(S)print(ans)