(真正的好東西)偏最小二乘回歸=多元線性回歸分析+典型相關(guān)分析+主成分分析復(fù)習(xí)課程

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

2、到許多限制。最典型的問題就是自變量之間的多重相關(guān)性。如果采用普通的最小二乘方法，這種變量多重相關(guān)性就會(huì)嚴(yán)重危害參數(shù)估計(jì)，擴(kuò)大模型誤差，并破壞模型的穩(wěn)定性。變量多重相關(guān)問題十分復(fù)雜，長(zhǎng)期以來(lái)在理論和方法上都未給出滿意的答案，這一直困擾著從事實(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)分

3、析+主成分分析由于偏最小二乘回歸在建模的同時(shí)實(shí)現(xiàn)了數(shù)據(jù)結(jié)構(gòu)的簡(jiǎn)化，因此，可以在二維平面圖上對(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è)因變量和p自變量。為了研究因變量和自變量的統(tǒng)計(jì)關(guān)系,我們觀測(cè)了n個(gè)樣本點(diǎn),由此構(gòu)成了自變量與因變量的數(shù)據(jù)

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

5、偏最小二乘回歸將通過(guò)實(shí)施對(duì)，的回歸,然后再表達(dá)成關(guān)于原變量，的回歸方程,k=1,2,q。1.2計(jì)算方法推導(dǎo)為了數(shù)學(xué)推導(dǎo)方便起見,首先將數(shù)據(jù)做標(biāo)準(zhǔn)化處理。X經(jīng)標(biāo)準(zhǔn)化處理后的數(shù)據(jù)矩陣記為=(，)，經(jīng)標(biāo)準(zhǔn)化處理后的數(shù)據(jù)矩陣記為=(，)。第一步記是的第一個(gè)成分，是的第一個(gè)軸，它是一個(gè)單位向量，既|=1。記是的第一個(gè)成分，=。是的第一個(gè)軸，并且|=1。如果要，能分別很好的代表X與Y中的數(shù)據(jù)變異信息，根據(jù)主成分分析原理，應(yīng)該有Var()maxVar()max另一方面，由于回歸建模的需要，又要求對(duì)有很大的解釋能力，有典型相關(guān)分析的思路，與的相關(guān)度應(yīng)達(dá)到最大值，既r（，）max因此，綜合起來(lái)，在偏最小二乘回歸

6、中，我們要求與的協(xié)方差達(dá)到最大，既Cov(，)=r(，)max正規(guī)的數(shù)學(xué)表述應(yīng)該是求解下列優(yōu)化問題，既s.t因此，將在|=1和|=1的約束條件下，去求()的最大值。如果采用拉格朗日算法，記s=(1)(1)對(duì)s分別求關(guān)于，和的偏導(dǎo)并令之為零，有=0(1-2)=0(1-3)=(1)=0(1-4)=(1)=0(1-5)由式(1-2)(1-5),可以推出記,所以,正是優(yōu)化問題的目標(biāo)函數(shù)值.把式(1-2)和式(1-3)寫成(1-6)(1-7)將式(1-7)代入式(1-6),有(1-8)同理,可得(1-9)可見,是矩陣的特征向量,對(duì)應(yīng)的特征值為.是目標(biāo)函數(shù)值,它要求取最大值,所以,是對(duì)應(yīng)于矩陣最大特征值的

7、單位特征向量.而另一方面,是對(duì)應(yīng)于矩陣最大特征值的單位特征向量.求得軸和后,即可得到成分然后,分別求和對(duì),的三個(gè)回歸方程(1-10)(1-11)(1-12)式中,回歸系數(shù)向量是(1-13)(1-14)(1-15)而,分別是三個(gè)回歸方程的殘差矩陣.第二步用殘差矩陣和取代和,然后,求第二個(gè)軸和以及第二個(gè)成分,有=是對(duì)應(yīng)于矩陣最大特征值的特征值,是對(duì)應(yīng)于矩陣最大特征值的特征向量.計(jì)算回歸系數(shù)因此,有回歸方程如此計(jì)算下去,如果的秩是,則會(huì)有（1-16）（1-17）由于,均可以表示成的線性組合,因此,式(1-17)還可以還原成關(guān)于的回歸方程形式，即k=1,2,q是殘差距陣的第k列。1.3交叉有效性下面要

8、討論的問題是在現(xiàn)有的數(shù)據(jù)表下,如何確定更好的回歸方程。在許多情形下,偏最小二乘回歸方程并不需要選用全部的成分進(jìn)行回歸建模,而是可以象在主成分分析一樣,采用截尾的方式選擇前m個(gè)成分,僅用這m個(gè)后續(xù)的成分就可以得到一個(gè)預(yù)測(cè)性較好的模型。事實(shí)上,如果后續(xù)的成分已經(jīng)不能為解釋提供更有意義的信息時(shí),采用過(guò)多的成分只會(huì)破壞對(duì)統(tǒng)計(jì)趨勢(shì)的認(rèn)識(shí),引導(dǎo)錯(cuò)誤的預(yù)測(cè)結(jié)論。在多元回歸分析一章中,我們?cè)谡{(diào)整復(fù)測(cè)定系數(shù)的內(nèi)容中討論過(guò)這一觀點(diǎn)。下面的問題是怎樣來(lái)確定所應(yīng)提取的成分個(gè)數(shù)。在多元回歸分析中,曾介紹過(guò)用抽樣測(cè)試法來(lái)確定回歸模型是否適于預(yù)測(cè)應(yīng)用。我們把手中的數(shù)據(jù)分成兩部分:第一部分用于建立回歸方程,求出回歸系數(shù)估計(jì)

9、量,擬合值以及殘差均方和;再用第二部分?jǐn)?shù)據(jù)作為實(shí)驗(yàn)點(diǎn),代入剛才所求得的回歸方程,由此求出。一般地,若有,則回歸方程會(huì)有更好的預(yù)測(cè)效果。若,則回歸方程不宜用于預(yù)測(cè)。在偏最小二乘回歸建模中,究竟應(yīng)該選取多少個(gè)成分為宜,這可通過(guò)考察增加一個(gè)新的成分后,能否對(duì)模型的預(yù)測(cè)功能有明顯的改進(jìn)來(lái)考慮。采用類似于抽樣測(cè)試法的工作方式,把所有n個(gè)樣本點(diǎn)分成兩部分:第一部分除去某個(gè)樣本點(diǎn)的所有樣本點(diǎn)集合(共含n-1個(gè)樣本點(diǎn)),用這部分樣本點(diǎn)并使用h個(gè)成分?jǐn)M合一個(gè)回歸方程;第二部分是把剛才被排除的樣本點(diǎn)代入前面擬合的回歸方程,得到在樣本點(diǎn)上的擬合值。對(duì)于每一個(gè)=1,2,n,重復(fù)上述測(cè)試,則可以定義的預(yù)測(cè)誤差平方和為,

10、有(1-18)定義Y的預(yù)測(cè)誤差平方和為,有(1-19)顯然,如果回歸方程的穩(wěn)健性不好,誤差就很大,它對(duì)樣本點(diǎn)的變動(dòng)就會(huì)十分敏感,這種擾動(dòng)誤差的作用,就會(huì)加大的值。另外,再采用所有的樣本點(diǎn),擬合含h個(gè)成分的回歸方程。這是,記第個(gè)樣本點(diǎn)的預(yù)測(cè)值為,則可以記的誤差平方和為,有(1-20)定義Y的誤差平方和為,有(1-21)一般說(shuō)來(lái),總是有大于,而則總是小于。下面比較和。是用全部樣本點(diǎn)擬合的具有h-1個(gè)成分的方程的擬合誤差;增加了一個(gè)成分,但卻含有樣本點(diǎn)的擾動(dòng)誤差。如果h個(gè)成分的回歸方程的含擾動(dòng)誤差能在一定程度上小于(h-1)個(gè)成分回歸方程的擬合誤差,則認(rèn)為增加一個(gè)成分,會(huì)使預(yù)測(cè)結(jié)果明顯提高。因此我們

11、希望的比值能越小越好。在SIMCA-P軟件中,指定即時(shí),增加成分就是有益的;或者反過(guò)來(lái)說(shuō),當(dāng)時(shí),就認(rèn)為增加新的成分,對(duì)減少方程的預(yù)測(cè)誤差無(wú)明顯的改善作用.另有一種等價(jià)的定義稱為交叉有效性。對(duì)每一個(gè)變量,定義(1-22)對(duì)于全部因變量Y,成分交叉有效性定義為(1-23)用交叉有效性測(cè)量成分對(duì)預(yù)測(cè)模型精度的邊際貢獻(xiàn)有如下兩個(gè)尺度。當(dāng)時(shí),成分的邊際貢獻(xiàn)是顯著的。顯而易見,與是完全等價(jià)的決策原則。對(duì)于k=1,2,q,至少有一個(gè)k,使得這時(shí)增加成分,至少使一個(gè)因變量的預(yù)測(cè)模型得到顯著的改善,因此,也可以考慮增加成分是明顯有益的。明確了偏最小二乘回歸方法的基本原理、方法及算法步驟后，我們將做實(shí)證分析。附錄

12、functionw=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);%functionc,m,v=norm1(C)%對(duì)數(shù)據(jù)進(jìn)行標(biāo)準(zhǔn)化處理n,s=size(C);fori=1:nforj=1:sc(i,j)=(C(i,j)-mean(C(:,j)/sqrt(cov(C(:,j);endendm=mean(C);forj=1:sv(1,j)=sqrt(cov(C(:,j);end%functiont,q,w,wh,f0,FF=fun717(px,py,C)%px自變

13、量的輸入個(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ù)y=c(:,px+1:px+py);%截取標(biāo)準(zhǔn)化的因變量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);%獲

14、得回歸系數(shù)p(:,1:px)=(E0*t(:,1)/(t(:,1)*t(:,1);fori=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);E(:,px*i+px+1:px*

15、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);endfors=1:pxq(:,s)=p(1,px*(s-1)+1:px*s);endn,d=size(q);forh=1:pxiw=eye(d);forj=1:h-1iw=iw*(eye(d)-w(:,j)*q(:,j);endwh(:,h)=iw*w(:,h);endforj=1:pyzr(j,:)=(regress1(y(:,j),t);%求回歸系數(shù)endforj=1:pxfori=1:py%生成標(biāo)準(zhǔn)化變量的方程的系數(shù)

16、矩陣w1=wh(:,1:j);zr1=(zr(i,1:j);f0(i,:,j)=(w1*zr1);endnormxy,meanxy,covxy=norm1(C);%normxy標(biāo)準(zhǔn)化后的數(shù)據(jù)矩陣%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);

17、FF(:,:,j)=fff,ff;%生成原始變量方程的常數(shù)項(xiàng)和系數(shù)矩陣end%functionr,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.2;RdYt=mean(Rdyt)form=1:pxRdYtt(1,m)=sum(RdYt(1,1:m);endforj=1:pyform=1:pyRdytt(j,m)=sum(

18、Rdyt(j,1:m);endendforj=1:pxform=1:pxRd(j,m)=RdYt(1,1:m)*(w(j,1:m).2);endendforj=1:pxVIP(j,:)=sqrt(px*ones(1,px)./RdYtt).*Rd(j,:);end%functionr,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

19、.2;RdXt=mean(Rdxt);form=1:pxRdXtt(1,m)=sum(RdXt(1,1:m);endforj=1:pxform=1:pxRdxtt(j,m)=sum(Rdxt(j,1:m);endend%forj=1:px%form=1:px%Rd(j,m)=RdXt(1,1:m)*(w(j,1:m).2);%end%end%forj=1:px%VIP(j,:)=sqrt(px*ones(1,px)./RdYtt).*Rd(j,:);%end%functiont,u=TU(px,py,C)%t提取的自變量的主成分%u提取的因變量的主成分c=norm1(C);y=c(:,px+1

20、: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);%functiondrew(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,HPZD,JPZD,TZ,ZG,GPK)gridonplot(YY(:,4),Y(:,4),+);lslin