第9講--多元線性回歸--主成分回歸

1、多元線性回歸主成分及其回歸第9講多元線性回歸解決的問(wèn)題系數(shù)矩陣系數(shù)矩陣 Y=XA建模：求解建模：求解回歸系數(shù)回歸系數(shù)A，該過(guò)程稱為建，該過(guò)程稱為建模模預(yù)報(bào)：在預(yù)報(bào)：在A已知時(shí)，對(duì)于新測(cè)已知時(shí)，對(duì)于新測(cè)Xnew，預(yù)報(bào)預(yù)報(bào)Ynew，稱為稱為預(yù)報(bào)預(yù)報(bào)例子某保健品含片產(chǎn)品，說(shuō)明書(shū)標(biāo)明：由營(yíng)養(yǎng)某保健品含片產(chǎn)品，說(shuō)明書(shū)標(biāo)明：由營(yíng)養(yǎng)物質(zhì)物質(zhì)A、B、C組成，產(chǎn)品標(biāo)注中寫出了每組成，產(chǎn)品標(biāo)注中寫出了每片中片中A、B、C物質(zhì)的含量。問(wèn)，如何認(rèn)定物質(zhì)的含量。問(wèn)，如何認(rèn)定？配置配置A、B、C的一組溶液，建立濃度與光的一組溶液，建立濃度與光吸收的關(guān)系。既建模求回歸系數(shù)吸收的關(guān)系。既建模求回歸系數(shù)將藥片配置成溶液，測(cè)吸光

2、，利用上面的將藥片配置成溶液，測(cè)吸光，利用上面的模型，預(yù)報(bào)濃度。模型，預(yù)報(bào)濃度。建模公式推導(dǎo)Y=XAXtY=XtXA(XtX)-1XtY=AE:學(xué)校教學(xué)pythonX.txtE:學(xué)校教學(xué)pythonY.txt問(wèn)題求解的關(guān)鍵步驟是什么？問(wèn)題求解的關(guān)鍵步驟是什么？方程數(shù)與未知數(shù)的關(guān)系設(shè)有規(guī)律上符合如下方程的一設(shè)有規(guī)律上符合如下方程的一 組實(shí)驗(yàn)數(shù)據(jù)組實(shí)驗(yàn)數(shù)據(jù) y= ax+b通過(guò)實(shí)驗(yàn)，不斷變更通過(guò)實(shí)驗(yàn)，不斷變更x，測(cè)得對(duì)應(yīng)的，測(cè)得對(duì)應(yīng)的y求求a,b的值，需要幾組這樣的數(shù)據(jù)？的值，需要幾組這樣的數(shù)據(jù)？唯一解唯一解最小二乘解最小二乘解y1y2ynx1 1x2 1 1xn 1ab=矩陣形式矩陣形式XtX是是

3、2*2的矩陣的矩陣方程數(shù)與未知數(shù)的關(guān)系設(shè)有規(guī)律上符合如下方程的一設(shè)有規(guī)律上符合如下方程的一 組實(shí)驗(yàn)數(shù)據(jù)組實(shí)驗(yàn)數(shù)據(jù) y= 1.2 x1 + 0.9 x2 + 3.3 x3通過(guò)實(shí)驗(yàn)，不斷變更通過(guò)實(shí)驗(yàn)，不斷變更x1、x2、x3，測(cè)得對(duì)應(yīng)的，測(cè)得對(duì)應(yīng)的y需要幾組這樣的數(shù)據(jù)？需要幾組這樣的數(shù)據(jù)？唯一解唯一解最小二乘解最小二乘解方程數(shù)小于未知數(shù)，一定無(wú)解嗎y= 1.2 x1 + 0.9 x2 + 3.3 x3當(dāng)當(dāng)X1,X2,X3存在線性相關(guān)時(shí)，問(wèn)題會(huì)怎樣？存在線性相關(guān)時(shí)，問(wèn)題會(huì)怎樣？如果如果x個(gè)數(shù)很多，樣本打不到要求，怎么辦個(gè)數(shù)很多，樣本打不到要求，怎么辦？現(xiàn)實(shí)中存在這樣的問(wèn)題嗎不同濃度成分相同的溶液，在

4、不同波長(zhǎng)不同濃度成分相同的溶液，在不同波長(zhǎng)x1、x2下下的吸光值的比值，溶液濃度變化，比值不變。的吸光值的比值，溶液濃度變化，比值不變。既既X1和和X2之間是線性相關(guān)的。之間是線性相關(guān)的。怎樣知道變量之間有相關(guān)性？答案：通過(guò)線性變化答案：通過(guò)線性變化主成分算法能解決這類問(wèn)題主成分算法能解決這類問(wèn)題死計(jì)算：檢查死計(jì)算：檢查XtX有沒(méi)有逆，沒(méi)逆，則線性相關(guān)有沒(méi)有逆，沒(méi)逆，則線性相關(guān)10主成份分析主成份分析 PCA Principle Component Analysis能有效的提取測(cè)量數(shù)據(jù)的有用信息能有效的提取測(cè)量數(shù)據(jù)的有用信息解決變量之間的相關(guān)性問(wèn)題解決變量之間的相關(guān)性問(wèn)題有效去除誤差，建立有效

5、的模型有效去除誤差，建立有效的模型11PCA分解算法原理分解算法原理采用非線性迭代偏最小二乘法采用非線性迭代偏最小二乘法(Nonlinear Iterative Partial Least Squares, NIPALS)方法方法分解量測(cè)矩陣分解量測(cè)矩陣S S = T Pt + E =tipi + E主成分示例主成分示例12方差最大方向方差最大方向NIPALS算法每次只求一個(gè)主成分，目前最大散差方向算法每次只求一個(gè)主成分，目前最大散差方向儀器的信噪比儀器的信噪比儀器測(cè)量時(shí)，信號(hào)強(qiáng)度要遠(yuǎn)遠(yuǎn)大于噪聲儀器測(cè)量時(shí)，信號(hào)強(qiáng)度要遠(yuǎn)遠(yuǎn)大于噪聲信號(hào)的數(shù)據(jù)的方差要遠(yuǎn)遠(yuǎn)大于噪聲的方差信號(hào)的數(shù)據(jù)的方差要遠(yuǎn)遠(yuǎn)大于噪聲

6、的方差所以，所以，PCA可以區(qū)別噪聲可以區(qū)別噪聲樣例樣例x0.91.10.80.87 22.21.92.1y1.21.00.92 1.11.81 1.91.72.5原數(shù)據(jù)圖PCA后15通過(guò)特征值比值判斷有效變量數(shù)通過(guò)特征值比值判斷有效變量數(shù)在在i/ i+i，應(yīng)該達(dá)到最大應(yīng)該達(dá)到最大值值根據(jù)根據(jù)i值，取值，取T和和P的前的前i列，即可扔掉噪聲列，即可扔掉噪聲16主成分回歸主成分回歸PCRPrinciple Component Regression是多元線性回歸！原來(lái) Y=XA 現(xiàn)在 Y=TAT為X的主成分得分，即X經(jīng)PCA分解后的得分因?yàn)門只是X的線性組合，提取了線性相關(guān)的部分，且只取前i列，所

7、以模型穩(wěn)定，去掉噪聲numpy中主成分中主成分分解分解SVD分解分解實(shí)矩陣的實(shí)矩陣的SVD(Singular Value Decomposition，奇，奇異值分解異值分解 )分解分解：分解結(jié)果：分解結(jié)果：A=USV其中其中S是對(duì)角矩陣是對(duì)角矩陣numpy中主成分中主成分分解分解-SVD程序代碼：程序代碼： B = np.linalg.svd(A,full_matrices=False)full_matrices=False一定要寫，否則會(huì)按復(fù)數(shù)一定要寫，否則會(huì)按復(fù)數(shù)分解分解分解結(jié)果：U=B0lamda=B1V = B2Lamda是所有的特征值，可以計(jì)算相鄰比值，決定主成分，它不是一個(gè)矩陣實(shí)例

8、光譜矩陣的SVD分解數(shù)據(jù)：E:學(xué)校教學(xué)教改項(xiàng)目教材數(shù)據(jù)S-093790.txt是一個(gè)16*6的矩陣看看能求解個(gè)特征值？16個(gè)？ 6個(gè)？96個(gè)？實(shí)例光譜矩陣的SVD分解data=np.mafromtxt(E:學(xué)校教學(xué)教改項(xiàng)目教材數(shù)據(jù)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)實(shí)例光譜矩陣的SVD

9、分解 ld = B1 for i in range(len(ld)-1):temp = ldi/ldi+1print (temp)4.964219451633.37725370576101.1942313561.470743841531.96976226926矩陣中有3個(gè)有效特征值根據(jù)有效特征值，設(shè)定PCA的得分和載荷實(shí)例光譜矩陣的SVD分解根據(jù)主成分，規(guī)劃得分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矩陣，特征值比值列表根據(jù)特征值比值，規(guī)劃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類調(diào)用p應(yīng)該先調(diào)用decompose方法，根據(jù)返回的特征之比值，確定主成分p再調(diào)用PCAdecompose方法，設(shè)定得分和載荷矩陣主成分回歸原來(lái) Y=XA現(xiàn)在 Y=TATtY=TtTATtT-1TtY=A求得最終的回歸系數(shù)：主成分 X=T

12、Pt因?yàn)镻是正交矩陣，所以 T=XPY=TA=XPA=XAnew主成分回歸數(shù)據(jù)E:學(xué)校教學(xué)pythonS-093843.txtE:學(xué)校教學(xué)pythonC-093843.txt求解方程 C=SAS是6*16的矩陣 所以StS的逆不存在S是光譜矩陣，光譜的不同波長(zhǎng)間線性相關(guān)，所以可以用PCR程序代碼建模過(guò)程S=np.mafromtxt(“E:學(xué)校教學(xué)pythonS-093843.txt )S=S.dataC=np.mafromtxt(“E:學(xué)校教學(xué)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程序代碼建模過(guò)程P=V.Tfor i in range(len(lmada)-1):temp = lamdai/lamdai+1print (temp)k=int(input(“主成分?jǐn)?shù)為:”)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)程序代碼預(yù)報(bào)對(duì)新測(cè)定Snew:C=np.dot(Snew, Al

14、ast)擴(kuò)展-能否用MLR、PCA類PCR類傳遞X，Y給PCRPCR內(nèi)，以X調(diào)用PCA，確定主成分?jǐn)?shù)根據(jù)確定的主成分?jǐn)?shù)，確定T、P，以T，Y建模，并結(jié)合P確定回歸系數(shù)建立預(yù)報(bào)方法。擴(kuò)展-能否用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擴(kuò)展-能否用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調(diào)用函數(shù)S=np.mafromt