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

2、問題十分復雜，長期以來在理論和方法上都未給出滿意的答案，這一直困擾著從事實際系統分析的工作人員。在偏最小二乘回歸中開辟了一種有效的技術途徑，它利用對系統中的數據信息進行分解和篩選的方式，提取對因變量的解釋性最強的綜合變量，辨識系統中的信息與噪聲，從而更好地克服變量多重相關性在系統建模中的不良作用。（3）偏最小二乘回歸之所以被稱為第二代回歸方法，還由于它可以實現多種數據分析方法的綜合應用。由于偏最小二乘回歸在建模的同時實現了數據結構的簡化，因此，可以在二維平面圖上對多維數據的特性進行觀察，這使得偏最小二乘回歸分析的圖形功能十分強大。在一次偏最小二乘回歸分析計算后，不但可以得到多因變量對多自變量的

3、回歸模型，而且可以在平面圖上直接觀察兩組變量之間的相關關系，以及觀察樣本點間的相似性結構。這種高維數據多個層面的可視見性，可以使數據系統的分析內容更加豐富，同時又可以對所建立的回歸模型給予許多更詳細深入的實際解釋。一、 偏最小二乘回歸的建模策略原理方法設有 q個因變量y1,.,yq和p自變量x1,.,xp。為了研究因變量和自變量的統計關系,我們觀測了n個樣本點,由此構成了自變量與因變量的數據表X=x1,.,xp和.Y=y1,.,yq。偏最小二乘回歸分別在X與Y中提取出成分t1 和u1 (也就是說, t1是x1,.,xp 的線形組合, u1是y1,.,yq 的線形組合).在提取這兩個成分時,為了

4、回歸分析的需要,有下列兩個要求:(1) t1和u1應盡可能大地攜帶他們各自數據表中的變異信息;(2) t1與u1 的相關程度能夠達到最大。這兩個要求表明，t1和u1 應盡可能好的代表數據表X和Y,同時自變量的成分t1對因變量的成分u1 又有最強的解釋能力。在第一個成分t1和 u1被提取后，偏最小二乘回歸分別實施X 對 t1的回歸以及 Y對u1 的回歸。如果回歸方程已經達到滿意的精度，則算法終止；否則,將利用 X被t1解釋后的殘余信息以及Y 被t1 解釋后的殘余信息進行第二輪的成分提取。如此往復，直到能達到一個較滿意的精度為止。若最終對 X共提取了 m個成分t1，tm，yk偏最小二乘回歸將通過實

5、施 yk對t1，tm， 的回歸,然后再表達成1.2計算方法推導 關于原變量x1 的回歸方程,k=1,2,q 。 ，xm，為了數學推導方便起見,首先將數據做標準化處理。X 經標準化處理后的數據矩陣記為E0=(E，E)np，Yj經標準化處理后的數據矩陣記為010pF0=(F01，F0q)np。00w1是E的第一個軸，第一步 記t1是E的第一個成分，它是一個單位向量，既|w1|=1。記u1是F0的第一個成分，u1=F0c。c 是F的第一個軸，并且|c|=1。 11012如果要t1，u1能分別很好的代表X與Y中的數據變異信息，根據主成分分析原理，應該有Var(u1)max Var(t1)max另一方面

6、，由于回歸建模的需要，又要求t1對u1有很大的解釋能力，有典型相關分析的思路，t1與u1的相關度應達到最大值，既r（t1，u1）max因此，綜合起來，在偏最小二乘回歸中，我們要求t1與u1的協方差達到最大，既Cov(t1，u1)=(t)Var(u)r(t111，u1) max正規的數學表述應該是求解下列優化問題，既因此，將在|w1|2=1和|c1|2=1的約束條件下，去求(w1E'0F值。如果采用拉格朗日算法，記s=w1E'0F0對s分別求關于w1，c1，1sw1sc1s1''c)的最大1c (w11'1w11)2 (c1c11)'和2的偏導并令

7、之為零，有=E'0F0c 211w1=0 (1 -2)=F'0E0w1'2c=0 (1-3)21=(w1w11)=0 (1-4)s2=(c1c11)=0 (1-5)'由式(1-2)(1-5),可以推出21=22=w1E0F0c1=<E0w1,F0c1>''記1=21=22=w'1E'0F0c1,所以,1正是優化問題的目標函數值. 把式(1-2)和式(1-3)寫成E0F0c1=1w1'(1-6) (1-7)F'E0w1=1c1將式(1-7)代入式(1-6),有''2E0F0F0E0w1=1

8、w1 (1-8)同理,可得F0E0E0F0c1=1c1 (1-9)''2可見,w1是矩陣E'0F0F'0E0的特征向量,對應的特征值為12.1是目標函數值,它要求取最大值,所以, w1是對應于E'0F0F'0E0矩陣最大特征值的單位特征向量.而另一方面, c1是對應于矩陣F'0E0E'0F0最大特征值1的單位特征向量. 求得軸w1和c1后,即可得到成分t1=E0w1u1=F0c12然后,分別求E0和F0對t1,u1的三個回歸方程式中,回歸系數向量是E0=t1p'1+E1(1-10) (1-11)F0=u1qF0=t1r&#

9、39;'1+F*11+F1(1-12)p1=E0t1|t1|'2'(1-13)q1=F0u1|u1|F0t1|t1|2'2(1-14)r1=(1-15)而E1,F*1,F1分別是三個回歸方程的殘差矩陣. 第二用殘差矩陣E1和F1取代E0和F0,然后,求第二個軸w2和c2以及第二個成分t2,u2,有t2=E1w2 u2=F1c2'2=<t2,u2>=w2E1F1c2'''2w是對應于矩陣E1F1F1E1最大特征值2的特征值, c2是對應于矩陣''p2=E1t2|t2|F1t2|t2|2'2'

10、;因此,有回歸方程 如此r2=E1=t2pF1=t2r''2+E2 +F22,如果X的秩是A,則會有E0=t1pF0=t1r''1+ +tAp+ +tAr''A（1-16）（1-17）1A+F由于,t1, ,tA均可以表示成E01, ,E0p的線性組合,因此,式(1-17)還可以還原成yk*=F0k關于xj*=E0k的回歸方程形式，即yk*=k1x1+ +kpxFAk*p+FAk k=1,2,q是殘差距陣FA的第k列。下面要討論的問題是在現有的數據表下,如何確定更好的回歸方程。在許多情形下,偏最小二乘回歸方程并不需要選用全部的成分t1, ,tA進

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

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

13、擬合值y的預測誤差平方和為PRESSPRESShj,有nhj=(yi=1ijhj(-i)2-y(1-18)定義Y 的預測誤差平方和為PRESSh,有pPRESSh=PRESSj=1hj(1-19)顯然,如果回歸方程的穩健性不好,誤差就很大,它對樣本點的變動就會十分敏感,這種擾動誤差的作用,就會加大PRESSh的值。另外,再采用所有的樣本點,擬合含h 個成分的回歸方程。這是,記第i個樣本hji,則可以記yj的誤差平方和為SShj,有 點的預測值為ynSShj=(yi=1ijhji)2-y(1-20)定義Y的誤差平方和為SSh,有pSSh=SS (1-21)hjj=1一般說來,總是有PRESSh大

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

15、預測誤差無明顯的改善作用.另有一種等價的定義稱為交叉有效性。對每一個變量yk,定義2Qhk=1-PRESShkSS(h-1)k(1-22)對于全部因變量Y,成分th交叉有效性定義為qQ=1-2hPRESSk=1hk=1-(h-1)kPRESSSS(h-1)h (1-23) SS用交叉有效性測量成分th對預測模型精度的邊際貢獻有如下兩個尺度。(1) 當Qh2(1-0.952)=0.0975時, th成分的邊際貢獻是顯著的。顯而易見, Qh20.0975與(PRESS(2) h/SSh-1)<0.952是完全等價的決策原則。 對于k=1,2,q,至少有一個k,使得這時增加成分th,至少使一個

16、因變量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(:,j); endendm=me

17、an(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); %norm1為標準化函數y=c(:,px+1:px+py); %截取標準化的因變量E0=c(:,1:px);F0=c(:,px+1:px+py);A=E0'*F0*

18、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(:,px*i+1:px*i+p

19、x);% 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(1,px*(s-1)+1:p

20、x*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); %normxy標準化后的數據矩陣

21、11 求回歸系數%meanxy每一列的均值%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)')-meanxy(1,px+1:px+py)'FF(:,:,j)=fff,ff; %生成原始變量方程的常數項和系數矩陣end% function r,Rdyt,RdYt,RdYtt,Rdytt,

22、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.2;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:pxRd(j,m)=RdYt(1

23、,1:m)*(w(j,1:m).2)'); 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.2;RdXt=mean(Rdxt);for m=1:pxRdXtt

24、(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).2)'); % 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 提取的因變量的主成分c=norm1(C);y=c(:,p

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