1、（真正的好東西）偏最小二乘回歸=多元線性回歸分析+典型相關分析+主成分分析偏最小二乘回歸是一種新型的 多元統計數據分析方法，它與1983年由伍德 和阿巴諾等人首次提出。近十年來，它在理論、方法和應用方面都得到了迅速 的發展。密西根大學的弗耐爾教授稱偏最小二乘回歸為第二代回歸分析方法。 偏最小二乘回歸方法在統計應用中的重要性主要的有以下幾個方面：（1）偏最小二乘回歸是一種 多因變量對多自變量 的回歸建模方法。（2）偏最小二乘回歸可以較好地解決許多以往用普通多元回歸無法解決的 問題。在普通多元線形回歸的應用中，我們常受到許多限制。最典型的問題就 是自變量之間的多重相關性。如果采用普通的最小二乘方法

2、，這種變量多重相 關性就會嚴重危害參數估計，擴大模型誤差，并破壞模型的穩定性。變量多重 相關問題十分復雜，長期以來在理論和方法上都未給出滿意的答案，這一直困 擾著從事實際系統分析的工作人員。在偏最小二乘回歸中開辟了一種有效的技 術途徑，它利用對系統中的數據信息進行分解和篩選的方式，提取對因變量的 解釋性最強的綜合變量，辨識系統中的信息與噪聲，從而更好地克服變量多重 相關性在系統建模中的不良作用。（3）偏最小二乘回歸之所以被稱為第二代回歸方法，還由于它可以實現多 種數據分析方法的綜合應用。偏最小二乘回歸二多元線性回歸分析+典型相關分析+主成分分析由于偏最小二乘回歸在建模的同時實現了數據結構的簡化

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

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

5、被ti解釋后的殘余信息以及丫被ti解釋后的殘余信息進行第二輪的成分 提取。如此往復，直到能達到一個較滿意的精度為止。若最終對 X 共提取了 m 個成分ti . tm偏最小二乘回歸將通過實施yk對ti. tm 的回歸撚后再表達成 y k 關于原變量 X i . X m 的回歸方程 ,k=i,2,. ,q 。1.2計算方法推導為了數學推導方便起見，首先將數據做標準化處理。X經標準化處理后的數 據矩陣記為Eo=（Eoi，Eop）n p，Yj經標準化處理后的數據矩陣記為F0 =（ F01，F0q ） n p 0第一步 記ti是Eo的第一個成分，Wi是Eo的第一個軸，它是一個單位向量，既 |W1| =

6、1o記U1是Fo的第一個成分，U1=FoC1 o c1是Fo的第一個軸，并且|c1|=1o如果要t1U1能分別很好的代表X與丫中的數據變異信息，根據主成分分析原理，應該有Var(uJ maxVar( t1)max另一方面，由于回歸建模的需要，又要求t1對U1有很大的解釋能力，有典型相 關分析的思路，t1與u1的相關度應達到最大值，既r （t1 u1）max因此，綜合起來，在偏最小二乘回歸中，我們要求 t1與U1的協方差達到最大,Cov(t1 uj=寸Var qjVar (小)r(t1uj max正規的數學表述應該是求解下列優化問題，既0 C1maxW1,C1s.tW 1 W 11C1 C11值

7、。如果采用拉格朗日算法，記s=W Eo Fo O -1(W w1 - 1)-2 (C1 C1 -1)對s分別求關于w1 C11和2的偏導并令之為零,由式(1-2)(1-5)，可以推出s '=Eo Fo C1 - 2 1 w1=ow1(1-2)S' lC=Fo Eo w1 - 22 C1=OCi(1-3)s=(w1 w1 1)=o1S'=-(C1 C1 1)=o22 12 2 w 1 E o FoC1Eow1 > FOC1(1-4)(1-5)w 1EoFoC1,所以，1正是優化問題的目標函數值.把式(1-2)和式(1-3)寫成IE o F0C1(1-6)F oEow

8、11C1(1-7)因此，將在|W1=1和|°=1的約束條件下，去求(W1 E。FocJ的最大(1-8)(1-9)將式(1-7)代入式(1-6),有E o Fo F o Eow1同理，可得F o EoE o Foc1!IQ可見,Wi是矩陣E 0F0F 0E0的特征向量，對應的特征值為1 . 1是目標函數值，它要 求取最大值，所以，wi是對應于E'oFoF '0E0矩陣最大特征值的單位特征向量而另 一方面,Ci是對應于矩陣F'oEoE'oFo最大特征值12的單位特征向量.求得軸Wi和Ci后,即可得到成分tiEoWiuiFoci然后,分別求Eo和Fo對ti,

9、 Ui的三個回歸方程EotiPiEi(i-iO)Fo1UiqiF i(i-ii)Fo1tir iFi(i-i2)式中，回歸系數向量是PiE oti(i-i3)lltiqi1F oui(i-i4)lluill2F oti(i-i5)rillti ll2而Ei,F i, Fi分別是三個回歸方程的殘差矩陣.第二步 用殘差矩陣Ei和Fi取代Eo和Fo，然后，求第二個軸W2和C2以及第 個成分t2,U2,有上2= Ei W2U2 = Fi C2112 t? ,u? w 2 E iFi C2W2是對應于矩陣E'lFiF'lEi最大特征值22的特征值,C2是對應于矩陣F'iEiE&#

10、39;iFi最大特征值的特征向量計算回歸系數P2E it2I|t2F it2因此，有回歸方程Elt2P2 E2如此計算下去，如果X的秩是A,則會有Eo tip'1tAp'A( 1-16)Fo "'1tA*AFa( 1-17)由于，如,tA均可以表示成Eo1,Eop的線性組合，因此，式(1-17)還可以還原成Fok關于Xj* Eok的回歸方程形式，即* *yk*k1X1kpX pFAkk=1,2,qFAk是殘差距陣FA的第k列。1.3交叉有效性下面要討論的問題是在現有的數據表下，如何確定更好的回歸方程。在許多情形下，偏最小二乘回歸方程并不需要選用全部的成分t1

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

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

13、,n,重復上述測試，則可以定義 yj 的預測誤差平方和為 PRESShj ,有nPRESShj(yijy?hj( i)2i1(1-18)定義 Y 的預測誤差平方和為PRESSh，有PRESShpPRESShjj1(1-19)顯然,如果回歸方程的穩健性不好 ,誤差就很大 ,它對樣本點的變動就會十分敏感這種擾動誤差的作用，就會加大PRESSh的值。另外,再采用所有的樣本點，擬合含h個成分的回歸方程。這是，記第i個樣本 點的預測值為?hji，則可以記yj的誤差平方和為sShj，有nSSij(yij ?iji )(1-20)i 1定義丫的誤差平方和為SSh，有pSShSShj(1-21)j i一般說來

14、,總是有PRESSh大于SSh，而SSh則總是小于SSh i。下面比較SSh i和PRESSh。 SSh i是用全部樣本點擬合的具有h-1個成分的方程的擬合誤差；PRESSh增加了一個成分th,但卻含有樣本點的擾動誤差。如果 h個成分的回歸方程的含擾動誤差能在一定程度上小于(h-1)個成分回歸方程的擬合誤差，則認為 增加一個成分th,會使預測結果明顯提高。因此我們希望(PRESSh/SSh J的比值 能越小越好。在SIMCA-P軟件中，指定(PRESSh /SSh 1)0.9521 '即.PRESS 0.95 SS 1時，增加成分th就是有益的；或者反過來說，當.PRESS 0.95.

15、 SSh 1時,就認為增加新的成分th,對減少方程的預測誤差無明顯 的改善作用.另有一種等價的定義稱為交叉有效性。對每一個變量yk，定義PRESShQhk 11)k(1-22)對于全部因變量丫,成分th交叉有效性定義為qPRESShk2 k 1hkPRESS,(1-23)Qh 1 1hSS；h 1)kSS(h1)用交叉有效性測量成分 th 對預測模型精度的邊際貢獻有如下兩個尺度。(1)當Qh (1 0.952) 0.0975時,th成分的邊際貢獻是顯著的。顯而易見,Q20.0975與(PRESSh/SShJ 0.952是完全等價的決策原則。 對于k=1,2，;q,至少有一個k,使得Qh20.0

16、975這時增加成分th,至少使一個因變量yk的預測模型得到顯著的改善，因此也 可以考慮增加成分th是明顯有益的。明確了偏最小二乘回歸方法的基本原理、方法及算法步驟后，我們將做實證分析。附錄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) %對數據進行標準化處理n,s=size(C);for i=1:nfor j=1:sc(i,j)=(C(i,j)-mean(C(:,j)/sqrt(cov(C

17、(:,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,C)% px自變量的輸入個數% py輸入因變量的個數。% C輸入的自變量和因變量組成的矩陣% t提取的主成分% q為回歸系數。% w最大特征值所對應的特征向量。% wh處理后的特征向量% f0回歸的標準化的方程系數% FF原始變量的回歸方程的系數c=norm1(C);%norm1y=c(:,px+1:px+py);%截取標準化的因變量E0=c(:,1:px);F0=c(:,px+1:px+py

18、);A=E0'*F0*F0'*E0;w(:,1)=maxdet(A);%求最大特征向量t(:,1)=E0*w(:,1);%提取主成分E(:,1:px)=E0-t(:,1)*(E0'*t(:,1)/(t(:,1)'*t(:,1)'%獲得回歸系數 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(:

19、,px*i+1:px*i+px);% maxdet 為求最大特征值的函數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)'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

20、(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)' %endfor j=1:pxfori=1:py %生成標準化變量的方程的系數矩陣w1=wh(:,1:j);zr1=(zr(i,1:j)'f0(i,:,j)=(w1*zr1)'endnormxy,meanxy,covxy=norm1(C);or

21、mxy 標準化后的數據矩陣%mea nx每一列的均值%covxy每一列的方差ccxx=ones(py,1)*meanxy(1,1:px);ccy=(covxy(1,px+1:px+py)'*ones(1,px);求回歸系數%nccx=ones(py,1)*(covxy(1,1:px);ff=ccy.*f0(:,:,j)./ccx;fff=-(sum(ccy.*ccxx.*f0(:,:,j)./ccx)')- meanxy(1,px+1:px+py)'FF(:,:,j)=fff,ff; % 生成 原始變量方程的常數項和系數矩陣end% function r,Rdyt,Rd

22、Yt,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:pyfor m=1:pyRdytt(j,m)=sum(Rdyt(j,1:m)');endendfor j=1:pxfor m=1

23、: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); r=r1(px+1:px+px,1:px)'Rdxt=r.A2;RdXt=mean(Rdx

24、t);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 提取的自變量的主成分% u 提取的因變量的主成分

25、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)'subplot(1,1,1,1)bar(f0(:,:,3)title('直方圖 ')legend('SG','TZBFB','FHL','JK',