(真正的好東西)偏最小二乘回歸=多元線性回歸分析+典型相關(guān)分析+主成分分析教學(xué)內(nèi)容

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

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

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

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

5、被ti解釋后的殘余信息以及Y被ti解釋后的殘余信息進(jìn)行第二輪的成分 提取。如此往復(fù)，直到能達(dá)到一個較滿意的精度為止。若最終對X 共提取了 m個成分ti tm偏最小二乘回歸將通過實施yk對ti tm 的回歸，然，后再表達(dá)成y k關(guān)于原變量x 1 x的回歸方程，k=i,2,q。X i X m，1.2計算方法推導(dǎo)為了數(shù)學(xué)推導(dǎo)方便起見，首先將數(shù)據(jù)做標(biāo)準(zhǔn)化處理。X經(jīng)標(biāo)準(zhǔn)化處理后的數(shù)據(jù)矩陣記為Eo=(Eo1 ,，E0P)np , Yj經(jīng)標(biāo)準(zhǔn)化處理后的數(shù)據(jù)矩陣記為F0 =( F01 ,，F(xiàn)0q ) n p。第一步 記ti是Eo的第一個成分，Wi是Eo的第一個軸，它是一個單位向量，既 l|w1H=1。記Ui是

6、Fo的第一個成分，U1 = F0 clo是Fo的第一個軸，并且|。|二1。如果要t1 U1能分別很好的代表X與Y中的數(shù)據(jù)變異信息，根據(jù)主成分分析原理，應(yīng)該有Var(U1) maxVar( t1) max另一方面，由于回歸建模的需要，又要求 G對5有很大的解釋能力，有典型相關(guān)分析的思路，t1與5的相關(guān)度應(yīng)達(dá)到最大值，既r (t1 U1) max因此，綜合起來，在偏最小二乘回歸中，我們要求 t1與U1的協(xié)方差達(dá)到最大，u1) maxCov(t1U1)= , Var(t)Var(5)r(t1正規(guī)的數(shù)學(xué)表述應(yīng)該是求解下列優(yōu)化問題，既maxW1，c1E o W1, Fo c1s.tW1 W11'

7、c1 c11值。因此，將在|wi|2二i和11cl |2二1的約束條件下，去求(w1 E0 F0 c1)的最大如果采用拉格朗日算法，記s=W Eo % C1 (Wi1 1)對s分別求關(guān)于W1c11和2的偏導(dǎo)并令之為零,2 (C1 C1-1)有wi=EoF0 C1 _ 2 1 W1=0(1-2)由式(1-2)(1-5)，可以推出s '7 = F0E0W1 _2 2c1=0 c1s , (W1 w1 - 1)=01S .'一一二一(C1 C1 -1)=02''2 12 2 w 1E 0F0c1E0W1,F0C1(1-3)(1-4)(1-5)記12 12 2 W1E0

8、F0C1,所以，1正是優(yōu)化問題的目標(biāo)函數(shù)值.把式(1-2)和式(1-3)寫成'E 0F0c11W1(1-6)一 F 0E0W11a(1-7)將式(1-7)代入式(1-6)，有2E 0 F0 F 0 E0W11 W1(1-8)同理，可得l' L l' L2F 0 E0E 0 F0G1 G(1-9). t. »» 9 、 可見,Wi是矩陣E 0F0F 0E0的特征向量，對應(yīng)的特征值為1 . 1是目標(biāo)函數(shù)值，它要求取最大值，所以，wi是對應(yīng)于E0FoF '0E0矩陣最大特征值的單位特征向量.而另一方面，Ci是對應(yīng)于矩陣F'oEoE'

9、oF。最大特征值12的單位特征向量.求得軸Wi和Ci后,即可得到成分tiEoWiuiFoci然后，分別求Eo和Fo對ti, Ui的三個回歸方程EotiPiEi(i-i0)Fouiqi(i-ii)FotiriFi(i-i2)式中，回歸系數(shù)向量是E otiPi211ti II2(i-i3)F oui qi 2lluill2(i-i4)F oti2|ti |2(i-i5)而Ei, F i,Fi分別是三個回歸方程的殘差矩陣.第二步 用殘差矩陣Ei和Fi取代Eo和Fo然后，求第二個軸W2和C2以及第個成分t2, U2,有t2= Ei W2U2 = Fi C2t2 ,u2w 2E iF1c2W2是對應(yīng)于矩

10、陣E'lFiF'lEi最大特征值2 2的特征值，C2是對應(yīng)于矩陣F'iEiE'iFi最大特征值的特征向量.計算回歸系數(shù),P2E it2211t2 if，F(xiàn) it2因此，有回歸方程El12 P 2 E2F1 t2r 2 F2如此計算下去，如果X的秩是A,則會有一 ''Eo tiP1 tAP A(1-16). '. '，一 一.Fo tir itAAFa(1-17)由于,ti, ,tA均可以表示成E01,Eop的線性組合，因此，式(1-17)還可以還原成ykFok關(guān)于Xj* Eok的回歸方程形式，即*yk*k1X1kpX pFAkk=

11、1,2,. ；qFAk是殘差距陣Fa的第k列。1.3交叉有效性下面要討論的問題是在現(xiàn)有的數(shù)據(jù)表下，如何確定更好的回歸方程。在許多情形下，偏最小二乘回歸方程并不需要選用全部的成分 t1, , tA進(jìn)行回3建模，而是可以象在主成分分析一樣，采用截尾的方式選擇前 m個成分(m A,A 秩(X),僅用這m個后續(xù)的成分t1, ,tm就可以得到一個預(yù)測性較好的模型。事實上，如 果后續(xù)的成分已經(jīng)不能為解釋F。提供更有意義的信息時，采用過多的成分只會破壞對統(tǒng)計趨勢的認(rèn)識,引導(dǎo)錯誤的預(yù)測結(jié)論。在多元回歸分析一章中,我們曾在調(diào)整復(fù)測定系數(shù)的內(nèi)容中討論過這一觀點。下面的問題是怎樣來確定所應(yīng)提取的成分個數(shù)。在多元回歸

12、分析中,曾介紹過用抽樣測試法來確定回歸模型是否適于預(yù)測應(yīng)用。我們把手中的數(shù)據(jù)分成兩部分:第一部分用于建立回歸方程,求出回歸系數(shù)估計量bB,擬合值兔以及殘差均方和?B ；再用第二部分?jǐn)?shù)據(jù)作為實驗點，代入剛才 所求得的回歸方程，由此求出外和?t o 一般地，若有?T ?B,則回歸方程會有更 好的預(yù)測效果。若?T2?B2 ,則回歸方程不宜用于預(yù)測。在偏最小二乘回歸建模中 ,究竟應(yīng)該選取多少個成分為宜,這可通過考察增加一個新的成分后,能否對模型的預(yù)測功能有明顯的改進(jìn)來考慮。采用類似于抽樣測試法的工作方式，把所有n個樣本點分成兩部分：第一部分除去某個樣本點i的所有樣本點集合(共含n-1個樣本點)，用這部

13、分樣本點并使用h個成分?jǐn)M合一個 回歸方程；第二部分是把剛才被排除的樣本點i代入前面擬合的回歸方程，得到y(tǒng)j 在樣本點i上的擬合值 商(”對于每一個i =1,2,n,重復(fù)上述測試，則可以定義 yj 的預(yù)測誤差平方和為PRESShj ,有n2PRESShj(yijy?hj( i)2(1-18)i1定義 Y 的預(yù)測誤差平方和為PRESSh ,有pPRESShPRESShj(1-19)j1顯然 ,如果回歸方程的穩(wěn)健性不好,誤差就很大,它對樣本點的變動就會十分敏感,這種擾動誤差的作用,就會加大PRESSh 的值。另外,再采用所有的樣本點，擬合含h個成分的回歸方程。這是，記第i個樣本點的預(yù)測值為 外，則可

14、以記yj的誤差平方和為SSj,有 nsshj(yj yhji )2(1-20)i 1定義Y的誤差平方和為SSh，有pSShSShj(1-21)j 1一般說來，總是有PRESSh大于SSh，而SSh則總是小于SSh 1。下面比較SSh 1和PRESSh。SSh 1是用全部樣本點擬合的具有h-1個成分的方程的擬合誤差；PRES0增加了一個成分th,但卻含有樣本點的擾動誤差。如果 h個成分的回歸方程的含擾動誤差能在一定程度上小于(h-1)個成分回歸方程的擬合誤差，則認(rèn)為增加一個成分th,會使預(yù)測結(jié)果明顯提高。因此我們希望(PRESS /SSh 1)的比值能越小越好。在SIMCA-P軟件中，指定(PR

15、ESSh / SSh 1)0.952即JPRESS 0.95南二時，增加成分th就是有益的；或者反過來說，當(dāng)JPRESS 0.95、SSh 1時,就認(rèn)為增加新的成分th,對減少方程的預(yù)測誤差無明顯 的改善作用.另有一種等價的定義稱為交叉有效性。對每一個變量 y定義Qhk 1PRES8SSh 1)k(1-22)PRESSh對于全部因變量Y,成分t h交叉有效性定義為QhqPRESShkk 1SS；h 1)kSS；h 1)(1-23)用交叉有效性測量成分th對預(yù)測模型精度的邊際貢獻(xiàn)有如下兩個尺度。(1) 當(dāng)Q2 (1 0.952) 0.0975時，th成分的邊際貢獻(xiàn)是顯著的。顯而易見，Q； 0.0

16、975與(PRESSh/S&i) 0.952是完全等價的決策原則。 對于k=1,2,q,至少有一個k,使得Qh20.0975這時增加成分th ,至少使一個因變量yk 的預(yù)測模型得到顯著的改善, 因此 ,也可以考慮增加成分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)%對數(shù)據(jù)進(jìn)行標(biāo)準(zhǔn)化處理n

17、,s=size(C);for i=1:nfor j=1:sc(i,j)=(C(i,j)-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,C)%px自變量的輸入個數(shù)%py輸入因變量的個數(shù)。%C輸入的自變量和因變量組成的矩陣%t提取的主成分%q為回歸系數(shù)。%w最大特征值所對應(yīng)的特征向量。%wh處理后的特征向量%f0回歸的標(biāo)準(zhǔn)化的方程系數(shù)%FF原始變量的回歸方程的系數(shù)c=norm1(C);%norm1為標(biāo)準(zhǔn)化函數(shù)截取

18、標(biāo)準(zhǔn)化的因變量求最大特征向量提取主成分%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(:,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)&#

19、39;*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)'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

20、(:,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:pxfori=1:py%生成標(biāo)準(zhǔn)化變量的方程的系數(shù)矩陣w1=wh(:,1:j);zr1=(zr(i,1:j)&

21、#39;f0(i,:,j)=(w1*zr1)'endnormxy,meanxy,covxy=norm1(C);%normxy 標(biāo)準(zhǔn)化后的數(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)')-meanxy(1,px+1:px+py)'FF(:,:

22、,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:pyfor m=1:pyRdytt(

23、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);r=

24、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% func

25、tion t,u=TU(px,py,C)%t 提取的自變量的主成分% 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)'subplot(1,1,1,1)bar(f0(:,:,3)title(' 直方圖 ')legend('SG','TZBFB','FHL','