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

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

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

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

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

5、被ti解釋后的殘余信息以及Y被ti解釋后的殘余信息進(jìn)行第二輪的成分 提取。如此往復(fù)，直到能達(dá)到一個(gè)較滿意的精度為止。若最終對(duì)X 共提取了 m個(gè)成分ti tm偏最小二乘回歸將通過實(shí)施yk對(duì)ti tm 的回歸，然，后再表達(dá)成y k關(guān)于原變量x 1 x的回歸方程，k=i,2,q。X i X m，1.2計(jì)算方法推導(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的第一個(gè)成分，Wi是Eo的第一個(gè)軸，它是一個(gè)單位向量，既 l|w1H=1。記Ui是

6、Fo的第一個(gè)成分，U1 = F0 clo是Fo的第一個(gè)軸，并且|。|二1。如果要t1 U1能分別很好的代表X與Y中的數(shù)據(jù)變異信息，根據(jù)主成分分析原理，應(yīng)該有Var(U1) maxVar( t1) max另一方面，由于回歸建模的需要，又要求 G對(duì)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)對(duì)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的特征向量，對(duì)應(yīng)的特征值為1 . 1是目標(biāo)函數(shù)值，它要求取最大值，所以，wi是對(duì)應(yīng)于E0FoF '0E0矩陣最大特征值的單位特征向量.而另一方面，Ci是對(duì)應(yīng)于矩陣F'oEoE'

9、oF。最大特征值12的單位特征向量.求得軸Wi和Ci后,即可得到成分tiEoWiuiFoci然后，分別求Eo和Fo對(duì)ti, Ui的三個(gè)回歸方程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分別是三個(gè)回歸方程的殘差矩陣.第二步 用殘差矩陣Ei和Fi取代Eo和Fo然后，求第二個(gè)軸W2和C2以及第個(gè)成分t2, U2,有t2= Ei W2U2 = Fi C2t2 ,u2w 2E iF1c2W2是對(duì)應(yīng)于矩

10、陣E'lFiF'lEi最大特征值2 2的特征值，C2是對(duì)應(yīng)于矩陣F'iEiE'iFi最大特征值的特征向量.計(jì)算回歸系數(shù),P2E it2211t2 if，F(xiàn) it2因此，有回歸方程El12 P 2 E2F1 t2r 2 F2如此計(jì)算下去，如果X的秩是A,則會(huì)有一 ''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個(gè)成分(m A,A 秩(X),僅用這m個(gè)后續(xù)的成分t1, ,tm就可以得到一個(gè)預(yù)測性較好的模型。事實(shí)上，如 果后續(xù)的成分已經(jīng)不能為解釋F。提供更有意義的信息時(shí)，采用過多的成分只會(huì)破壞對(duì)統(tǒng)計(jì)趨勢的認(rèn)識(shí),引導(dǎo)錯(cuò)誤的預(yù)測結(jié)論。在多元回歸分析一章中,我們?cè)谡{(diào)整復(fù)測定系數(shù)的內(nèi)容中討論過這一觀點(diǎn)。下面的問題是怎樣來確定所應(yīng)提取的成分個(gè)數(shù)。在多元回歸

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

13、分樣本點(diǎn)并使用h個(gè)成分?jǐn)M合一個(gè) 回歸方程；第二部分是把剛才被排除的樣本點(diǎn)i代入前面擬合的回歸方程，得到y(tǒng)j 在樣本點(diǎn)i上的擬合值 商(”對(duì)于每一個(gè)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)健性不好,誤差就很大,它對(duì)樣本點(diǎn)的變動(dòng)就會(huì)十分敏感,這種擾動(dòng)誤差的作用,就會(huì)加大PRESSh 的值。另外,再采用所有的樣本點(diǎn)，擬合含h個(gè)成分的回歸方程。這是，記第i個(gè)樣本點(diǎn)的預(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是用全部樣本點(diǎn)擬合的具有h-1個(gè)成分的方程的擬合誤差；PRES0增加了一個(gè)成分th,但卻含有樣本點(diǎn)的擾動(dòng)誤差。如果 h個(gè)成分的回歸方程的含擾動(dòng)誤差能在一定程度上小于(h-1)個(gè)成分回歸方程的擬合誤差，則認(rèn)為增加一個(gè)成分th,會(huì)使預(yù)測結(jié)果明顯提高。因此我們希望(PRESS /SSh 1)的比值能越小越好。在SIMCA-P軟件中，指定(PR

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

16、975與(PRESSh/S&i) 0.952是完全等價(jià)的決策原則。 對(duì)于k=1,2,q,至少有一個(gè)k,使得Qh20.0975這時(shí)增加成分th ,至少使一個(gè)因變量yk 的預(yù)測模型得到顯著的改善, 因此 ,也可以考慮增加成分th 是明顯有益的。明確了偏最小二乘回歸方法的基本原理、方法及算法步驟后，我們將做實(shí)證分析。附錄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)%對(duì)數(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自變量的輸入個(gè)數(shù)%py輸入因變量的個(gè)數(shù)。%C輸入的自變量和因變量組成的矩陣%t提取的主成分%q為回歸系數(shù)。%w最大特征值所對(duì)應(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ù)項(xiàng)和系數(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','