(真正地好東西)偏最小二乘回歸=多元線性回歸分析報告+典型相關分析報告+主成分分析報告

1、實用標準文檔偏最小二乘回歸是一種新型的 多元統(tǒng)計數(shù)據(jù)分析方法，它與1983年由伍德 和阿巴諾等人首次提出。近十年來，它在理論、方法和應用方面都得到了迅速的 發(fā)展。密西根大學的弗耐爾教授稱偏最小二乘回歸為第二代回歸分析方法。 偏最小二乘回歸方法在統(tǒng)計應用中的重要性主要的有以下幾個方面：(1)偏最小二乘回歸是一種 多因變量對多自變量 的回歸建模方法。(2)偏最小二乘回歸可以較好地解決許多以往用普通多元回歸無法解決的 問題。在普通多元線形回歸的應用中， 我們常受到許多限制。最典型的問題就是 自變量之間的多重相關性。如果采用普通的最小二乘方法，這種變量多重相關性 就會嚴重危害參數(shù)估計，擴大模型誤差，并

2、破壞模型的穩(wěn)定性。變量多重相關問 題十分復雜，長期以來在理論和方法上都未給出滿意的答案， 這一直困擾著從事 實際系統(tǒng)分析的工作人員。在偏最小二乘回歸中開辟了一種有效的技術途徑，它利用對系統(tǒng)中的數(shù)據(jù)信息進行分解和篩選的方式， 提取對因變量的解釋性最強的 綜合變量，辨識系統(tǒng)中的信息與噪聲，從而更好地克服變量多重相關性在系統(tǒng)建 模中的不良作用。(3)偏最小二乘回歸之所以被稱為第二代回歸方法，還由于它可以實現(xiàn)多 種數(shù)據(jù)分析方法的綜合應用。偏最小二乘回歸二多元線性回歸分析+典型相關分析+主成分分析由于偏最小二乘回歸在建模的同時實現(xiàn)了數(shù)據(jù)結(jié)構的簡化，因此，可以在二維平面圖上對多維數(shù)據(jù)的特性進行觀察，這使得

3、偏最小二乘回歸分析的圖形功能 十分強大。在一次偏最小二乘回歸分析計算后，不但可以得到多因變量對多自變 量的回歸模型，而且可以在平面圖上直接觀察兩組變量之間的相關關系，以及觀察樣本點間的相似性結(jié)構。這種高維數(shù)據(jù)多個層面的可視見性，可以使數(shù)據(jù)系統(tǒng) 的分析內(nèi)容更加豐富，同時又可以對所建立的回歸模型給予許多更詳細深入的實 際解釋。偏最小二乘回歸的建模策略原理方法1.1 建模原理設有q個因變量 yi,.,yq和p自變量 xh.,xp。為了研究因變量和自變 量的統(tǒng)計關系，我們觀測了n個樣本點，由此構成了自變量與因變量的數(shù)據(jù)表X=xi,.,xp和.Y= yi,.,yq。偏最小二乘回歸分別在X與Y中提取出成分

4、ti和 u1(也就是說，t1是x1,.,xp的線形組合，u1是y1,.,yq的線形組合).在提取這兩個成分時，為了回歸分析的需要，有下列兩個要求：(1) ti和Ui應盡可能大地攜帶他們各自數(shù)據(jù)表中的變異信息；(2) ti與ui的相關程度能夠達到最大。這兩個要求表明，ti和ui應盡可能好的代表數(shù)據(jù)表X和Y,同時自變量的成分ti對因變量的成分Ui又有最強的解釋能力。在第一個成分ti和ui被提取后，偏最小二乘回歸分別實施 X對ti的回歸以 及Y對ui的回歸。如果回歸方程已經(jīng)達到滿意的精度，則算法終止；否則 ，將 利用X被ti解釋后的殘余信息以及Y被ti解釋后的殘余信息進行第二輪的成分 提取。如此往復

5、，直到能達到一個較滿意的精度為止。若最終對X共提取了 m個成分titm 偏最小二乘回歸將通過實施yk對titm的回歸，然后再表達成yk關于原變量x i . Xm的回歸方程，k=i,2,q。1.2 計算方法推導為了數(shù)學推導方便起見，首先將數(shù)據(jù)做標準化處理。X經(jīng)標準化處理后的數(shù) 據(jù)矩陣記為Eo=(Ei,，Eop)n刈，Yj經(jīng)標準化處理后的數(shù)據(jù)矩陣記為F0 =( F0i ,，F(xiàn)0q ) nxp 第一步 記ti是Eo的第一個成分，Wi是Eo的第一個軸，它是一個單位向量,既| wi|=i 。記Ui是Fo的第一個成分，u產(chǎn)Fo。是Fo的第一個軸，并且| Ci|=i如果要ti Ui能分別很好的代表X與Y中的

6、數(shù)據(jù)變異信息，根據(jù)主成分分析原理，應該有Var( u1)maxVar( ti) max另一方面，由于回歸建模的需要，又要求ti對Ui有很大的解釋能力，有典型相關 分析的思路，3與5的相關度應達到最大值，既r( tiUi) t max因此，綜合起來，在偏最小二乘回歸中，我們要求 ti與Ui的協(xié)方差達到最大，既Cov(tiui)= JVar (tJVar (uJ r( tiui) t max正規(guī)的數(shù)學表述應該是求解下列優(yōu)化問題，既W去求(W Eo FO)Ci。一 i)s.t因此，將在| w/l 2=i和| Ci| 2二i的約束條件下,的最大值。如果采用拉格朗日算法，記s=W Eo Fo Ci %

7、i( Wi Wi i)九 2(對s分別求關于WiCi九1和九2的偏導并令之為零,：s _而一EoFo C12 i wi=0(i -2)s _ F Eo wi .ciCi=0(i-3)文案大全(i-4)-s ,1才一(Wi Wi一1)二0-i*=-(CCiT)=0(1-5)由式(1-2)(1-5),可以推出2.i 2 .i, 2 二 w 1E 0 F0 Ci = ； E 0 v1，F(xiàn)0C1記& =2% =2九2 =w1E0F0C1,所以，d正是優(yōu)化問題的目標函數(shù)值.把式(1-2)和式(1-3)寫成 E 0F0c1 - T1W1(1-6) F 0 E 0W1 = 口心(1-7)將式(1-7)代入式

8、(1-6),有 L2E 0 F0 F 0 E0w1 =11 w1(1-8)同理,可得l L l L. 2F 0E0E 0 F0c1 - u1 C1(1-9)可見，W1是矩陣E0F0F0E0的特征向量，對應的特征值為日：區(qū)是目標函數(shù)值，它要求取最大值,所以，W1是對應于E0F0F0E0矩陣最大特征值的單位特征向量.而另一方面，C1是對應于矩陣F0E0E0F0最大特征值32的單位特征向量.求得軸W1和g后，即可得到成分t - E0W1U1 = F0C1然后，分別求E。和F。對G, 5的三個回歸方程E。=Lp1E1(1-10)F0 一四一 F 1(1-11) , . F0 =tj 1F1(1-12)

9、式中，回歸系數(shù)向量是pl 二,E oti而q1 二F oUi2l|Ui II2ri一F ot12lltill2(1-13)(1-14)(1-15)而E1, F*1, F1分別是三個回歸方程的殘差矩陣.第二步 用殘差矩陣E1和R取代Eo和Fo,然后，求第二個軸W2和C2以及第二個成分t2, u2,有t2 = E1 w2U2 = F1 C22 二:二12, U2 二 w2 E1F1C2. 、 一 9 , 、W2是對位于矩陣 E1F1F 1E1最大特征值 外的特征值,C2是對位于矩陣 . . 、. .F1E1E1F1最大特征值的特征向量.計算回歸系數(shù),E 1t2p2 :211t2 一 F 1t2 r

10、2211t2因此，有回歸方程E1 = 12 p 2 E2F1 = 12r 2 F2如此計算下去，如果X的秩是A,則會有.一一.Eo=t1p1+tApA(1-16) . . . Fo =t-+ +tAa+Fa(1-17)由于，t1,，tA均可以表示成Eo1,，Eop的線性組合，因此，式(1-17)還可以還原成yk =Fk關于Xj* = Ek的回歸方程形式，即*y=c(kiXi + +kpXp + FAk k=1,2,，qFAk是殘差距陣Fa的第k歹I。1.3交叉有效性下面要討論的問題是在現(xiàn)有的數(shù)據(jù)表下，如何確定更好的回歸方程。在許多 情形下，偏最小二乘回歸方程并不需要選用全部的成分tl,，tA進

11、行回歸建模，而是可以象在主成分分析一樣，采用截尾的方式選擇前m個成分 (m A,A=秩(X),僅用這m個后續(xù)的成分ti,，tm就可以得到一個預測性較好 的模型。事實上，如果后續(xù)的成分已經(jīng)不能為解釋 Fo提供更有意義的信息時，采 用過多的成分只會破壞對統(tǒng)計趨勢的認識，引導錯誤的預測結(jié)論。在多元回歸分 析一章中，我們曾在調(diào)整復測定系數(shù)的內(nèi)容中討論過這一觀點。下面的問題是怎樣來確定所應提取的成分個數(shù)。在多元回歸分析中，曾介紹過用抽樣測試法來確定回歸模型是否適于預測應 用。我們把手中的數(shù)據(jù)分成兩部分：第一部分用于建立回歸方程，求出回歸系數(shù)估 計量bB ,擬合值%以及殘差均方和?B ;再用第二部分數(shù)據(jù)作

12、為實驗點，代入剛才 所求得的回歸方程，由此求出？t和赍。一般地,若有夕；歸田；,則回歸方程會有更 好的預測效果。若 有次B ,則回歸方程不宜用于預測。在偏最小二乘回歸建模中，究竟應該選取多少個成分為宜，這可通過考察增 加一個新的成分后，能否對模型的預測功能有明顯的改進來考慮。采用類似于抽 樣測試法的工作方式，把所有n個樣本點分成兩部分：第一部分除去某個樣本點i 的所有樣本點集合(共含n-1個樣本點)，用這部分樣本點并使用h個成分擬合一 個回歸方程；第二部分是把剛才被排除的樣本點i代入前面擬合的回歸方程，得到 yj在樣本點i上的擬合值?可)。對于每一個i=1,2,n,重復上述測試，則可以 定義y

13、j的預測誤差平方和為PRESShj,有nPRESShj = (yj ?hj(，)2(1-18)i 1定義Y的預測誤差平方和為PRESSh,有pPRESS =Z PRESSj(1-19)j=i顯然，如果回歸方程的穩(wěn)健性不好，誤差就很大，它對樣本點的變動就會十分敏感 這種擾動誤差的作用,就會加大PRESSh的值。另外，再采用所有的樣本點，擬合含h個成分的回歸方程。這是，記第i個樣 本點的預測值為？加，則可以記yj的誤差平方和為SShj,有(1-20)nSGj(yj -篩)2i 4定義Y的誤差平方和為SSh,有pSS 八 SSj(1-21)j 1般說來，總是有PRES與大于SSh,而SSh則總是小于

14、SSh_1。下面比較SSh和PRESSh。S&是用全部樣本點擬合的具有h-1個成分的方程的擬合誤差PRESSh增加了一個成分th，但卻含有樣本點的擾動誤差。如果h個成分的回歸方程的含擾動誤差能在一定程度上小于 (h-1)個成分回歸方程的擬合誤差，則認為 增加一個成分th,會使預測結(jié)果明顯提高。因此我們希望 (PRESSh / SSh)的比值 能越小越好。在SIMCA-P軟件中，指定(PRESSh / SSh) 0.957SSh二時，就認為增加新的成分th,對減少方程的預測誤差無明 顯的改善作用.另有一種等價的定義稱為交叉有效性。對每一個變量y一定義2 hkd PRESS-1 SSh 4)k(1

15、-22)對于全部因變量Y,成分th交叉有效性定義為“ PRESSk 1S SS(h1)k1 PRESShSSh)(1-23)用交叉有效性測量成分th對預測模型精度的邊際貢獻有如下兩個尺度。(1) 當Q：(1 -0.952)=0.0975時，th成分的邊際貢獻是顯著的。顯而易見，Q；*0.0975與(PRESSh/S&)0.952是完全等價的決策原貝限 對于k=1,2,q,至少有一個k,使得一 2 一一一 Qh -0.0975這時增加成分th,至少使一個因變量yk的預測模型得到顯著的改善，因此, 也可以考慮增加成分量是明顯有益的。明確了偏最小二乘回歸方法的基本原理、方法及算法步驟后，我們將做實證

16、分析。function w=maxdet(A)%求矩陣的最大特征值v,d=eig(A);n,p=size(d);d1=d*ones(p,1);d2=max(d1);i=find(d1=d2);w=v(：,i)；%function c,m,v=norm1(C)%寸數(shù)據(jù)進行標準化處理n,s=size(C);for i=1:nfor j=1:sc(i,j)=(C(iJ)-mean(C(：,j)/sqrt(cov(C(：,j)；endendm=mean(C);for j=1:sv(1,j)=sqrt(cov(C(:,j);end%function t,q,w,wh,f0,FF=fun717(px,py

17、,C)%px自變量的輸入個數(shù)%py輸入因變量的個數(shù)。%C輸入的自變量和因變量組成的矩陣% t提取的主成分% q為回歸系數(shù)。% w最大特征值所對應的特征向量。% wh處理后的特征向量% f0回歸的標準化的方程系數(shù)% FF原始變量的回歸方程的系數(shù)c=norm1(C);%norm1為標準化函數(shù)y=c(:,px+1:px+py);%截取標準化的因變量E0=c(:,1:px);F0=c(:,px+1:px+py);A=E0*F0*F0*E0;w(:,1)=maxdet(A);%t(:,1)=E0*w(:,1);%求最大特征向量提取主成分E(:,1:px尸E0-t(:,1)*(E0*t(:,1)/(t(:

18、,1)*t(:,1);% 獲得回歸系數(shù)p(:,1:px)=(E0*t(:,1)/(t(:,1)*t(:,1);for i=0:px-2B(:,px*i+1:px*i+px)=E(:,px*i+1:px*i+px)*F0*F0*E(:,px*i+1:px*i+px) ；w(:,i+2)=maxdet(B(:,px*i+1:px*i+px);% maxdet為求最大特征值的函數(shù)t(:,i+2)=E(:,px*i+1:px*i+px)*w(:,i+2);p(:,px*i+px+1:px*i+2*px)=(E(:,px*i+1:px*i+px)*t(:,i+2)/(t(:,i+2)*t(:,i+2)；

19、E(:,px*i+px+1:px*i+2*px)=E(:,px*i+1:px*i+px)-t(:,i+2)*(E(:,px*i+1:px* i+px)*t(:,i+2)/(t(:,i+2)*t(:,i+2);endfor s=1:pxq(:,s)=p(1,px*(s-1)+1:px*s);endn,d=size(q);for h=1:pxiw=eye(d);for j=1:h-1iw=iw*(eye(d)-w(:,j)*q(:,j);endwh(:,h)=iw*w(:,h);endfor j=1:pyzr(j,:)=(regress1(y(:,j),t);%求回歸系數(shù)endfor j=1:px

20、fori=1:py%生成標準化變量的方程的系數(shù)矩陣w1=wh(:,1:j);zr1=(zr(i,1:j);f0(i,:,j)=(w1*zr1);endnormxy,meanxy,covxy=norm1(C);%normxy標準化后的數(shù)據(jù)矩陣%meanxyj一列的均值%covxy每一列的方差ccxx=ones(py,1)*meanxy(1,1:px);ccy=(covxy(1,px+1:px+py)*ones(1,px);ccx=ones(py,1)*(covxy(1,1:px);ff=ccy*f0(:,:,j)./ccx;fff=-(sum(ccy.*ccxx.*f0(:,:,j)./ccx)

21、-meanxy(1,px+1:px+py);FF(:,:,j)=fff,ff;%生成原始變量方程的常數(shù)項和系數(shù)矩陣end%function r,Rdyt,RdYt,RdYtt,Rdytt,VIP=fun8y(px,py,c)X=c(:,1:px);Y=c(:,px+1:px+py);x=norm1(X);y=norm1(Y);t,q,w=fun717(px,py,X,Y);r1=corrcoef(y,t);r=r1(py+1:px+py,1:py);Rdyt=r.A2;RdYt=mean(Rdyt)for m=1:pxRdYtt(1,m)=sum(RdYt(1,1:m);endfor j=1：

22、pyfor m=1:pyRdytt(j,m)=sum(Rdyt(j,1:m);endendfor j=1:pxfor m=1:pxRd(j,m尸RdYt(1,1:m)*(w(j,1:m)A2);endendfor j=1:pxVIP(j,:)=sqrt(px*ones(1,px)./RdYtt)*Rd(j,:);end%function r,Rdxt,RdXt,RdXtt,Rdxtt=fun8x(px,py,c)X=c(:,1:px);Y=c(:,px+1:px+py);x=norm1(X);y=norm1(Y);t,q,w=fun717(px,py,X,Y);r1=corrcoef(x,t)

23、;r=r1(px+1:px+px,1:px);Rdxt=r.A2;RdXt=mean(Rdxt);for m=1:pxRdXtt(1,m)=sum(RdXt(1,1:m);endfor j=1:pxfor m=1:pxRdxtt(j,m)=sum(Rdxt(j,1:m);endend% for j=1:px% for m=1:px%Rd(j,m尸RdXt(1,1:m)*(w(j,1:m)A2);% end% end% for j=1:px% VIP(j,:)=sqrt(px*ones(1,px)./RdYtt)*Rd(j,:);% end%function t,u=TU(px,py,C)%t提

24、取的自變量的主成分%u提取的因變量的主成分c=norm1(C);y=c(:,px+1:px+py);E0=c(:,1:px);F0=c(:,px+1:px+py);A=E0*F0*F0*E0;w(:,1)=maxdet(A);t(:,1)=E0*w(:,1);B=F0*E0*E0*F0;cc(:,1)=maxdet(B);u(:,1)=F0*cc(:,1);%function drew(px,py,c)X=c(:,1:px);Y=c(:,px+1:px+py);line,l=size(Y);t,q,w,wh,f0,FF=fun717(px,py,c);YY=X*FF(:,2:px+1,3)+ones(line,1)*FF(:,1,3