算法大全第30章__偏最小二乘回歸_h

1、 第三十章偏最小二乘回歸在實際問題中，經常遇到需要研究兩組多重相關變量間的相互依賴關系，并研究用一組變量（常稱為自變量或預測變量）去預測另一組變量（常稱為因變量或響應變量）， 除了最小二乘準則下的經典多元線性回歸分析（MLR），提取自變量組主成分的主成分回歸分析（PCR）等方法外，還有近年發(fā)展起來的偏最小二乘（PLS）回歸方法。 偏最小二乘回歸提供一種多對多線性回歸建模的方法，特別當兩組變量的個數(shù)很多，且都存在多重相關性，而觀測數(shù)據(jù)的數(shù)量（樣本量）又較少時，用偏最小二乘回歸建立的模型具有傳統(tǒng)的經典回歸分析等方法所沒有的優(yōu)點。 偏最小二乘回歸分析在建模過程中集中了主成分分析，典型相關分析和線性回

2、歸分析方法的特點，因此在分析結果中，除了可以提供一個更為合理的回歸模型外，還可以同時完成一些類似于主成分分析和典型相關分析的研究內容，提供更豐富、深入的一些信息。 本章介紹偏最小二乘回歸分析的建模方法；通過例子從預測角度對所建立的回歸模型進行比較。 1 偏最小二乘回歸考慮 p 個變量 y1 , y2 , yp 與 m 個自變量 x1 , x2 , xm 的建模問題。偏最小二乘 回歸的基本作法是首先在自變量集中提出第一成分t1 （ t1 是 x1, xm 的線性組合， 且 盡可能多地提取原自變量集中的變異信息）；同時在因變量集中也提取第一成分u1 ，并要求t1 與u1 相關程度達到最大。然后建立

3、因變量 y1, yp 與t1 的回歸，如果回歸方 程已達到滿意的精度，則算法中止。否則繼續(xù)第二對成分的提取，直到能達到滿意的精 度為止。若最終對自變量集提取 r 個成分t1 , t2 , tr ，偏最小二乘回歸將通過建立 y1, yp 與t1,t2,tr 的回歸式，然后再表示為 y1, yp 與原自變量的回歸方程式， 為了方便起見，不妨假定 p 個因變量 y1 , yp 與 m 個自變量 x1 , xm 均為標準 化變量。因變量組和自變量組的n 次標準化觀測數(shù)據(jù)陣分別記為 x1m y11y1p x11= = FE#， ynp #yn100xn1xnm 偏最小二乘回歸分析建模的具體步驟如下：-5

4、31- （1）分別提取兩變量組的第一對成分，并使之相關性達最大。假設從兩組變量分別提出第一對成分為t1 和u1，t1 是自變量集 X = (x1, x m)T的 線性組合：t1 = w11x1 + + w1m xm = w X ，u 是因變量集Y = ( y , y ) 的線性組 TT111p合： u1= v 1y1 +1 + v y 1=p vp TY 。1為了回歸分析的需要，要求： t1 和u1 各自盡可能多地提取所在變量組的變異信息； t1 和u1 的相關程度達到最大。 由兩組變量集的標準化觀測數(shù)據(jù)陣 E0 和 F0 ，可以計算第一對成分的得分向量，記 為t1 和u1 ：x11x1m w

5、11t11 = = #t = E w#10 1xn1xnm w1m tn1y11y1 pv11 u11 = =u = F v#yn110 1ynp v1pun1第一對成分t1 和u1 的協(xié)方差Cov(t1, u1) 可用第一對成分的得分向量t1 和u1 的內積來計算。故而以上兩個要求可化為數(shù)學上的條件極值問題： = wT ET F x max110 10 110 0 122Tw = 1,vv = 1ww1Tv111 1v ，使q = w E F v 最大。問 TT利用Lagrange乘數(shù)法，問題化為求單位向量 w 和11110 0 1題的求解只須通過計算m m 矩陣 M = ET F FT E

6、 的特征值和特征向量，且 M 的最大特0 0 001征值為q ，相應的單位特征向量就是所求的解 2w 計算得到v =w ，而v 可由 FTE w 。q1111100 11（2）建立 y1 , yp 對t1 的回歸及 x1 , xm 對t1 的回歸。假定回歸模型為 -532- E = t a + ET01 11F = u b+ FT01 11其中a1 = (a11,a1m ) ， b = (b , b ) 分別是多對一的回歸模型中的參數(shù)向量， TT1111 pE1 和 F1 是殘差陣。回歸系數(shù)向量a1, b1 的最小二乘估計為 2= ET tt10 11，2b1 = FTtt0 11稱a1, b

7、1為模型效應負荷量。（3）用殘差陣 E1 和 F1代替 E0 和 F0 重復以上步驟。 記 E= t a ，F(xiàn) = t b，則殘差陣 ，F(xiàn) = F - F。如果殘差陣 FE = E - ETT01 101 11001001中元素的絕對值近似為 0，則認為用第一個成分建立的回歸式精度已滿足需要了，可以 停止抽取成分。否則用殘差陣 E1 和 F1 代替 E0 和 F0 重復以上步驟即得： w = (w； v= (v , v分 別為第 二 對成分 的 權數(shù)。 而 )T)T, w2212m2212 pt2 = E1w2 ， u2 = F1v2 為第二對成分的得分向量。ET tFT t1 2 a =b

8、=分別為 X ,Y 的第二對成分的負荷量。這時有 1 2，2222tt22E = t a + t a + ETT01 12 22F = tb+ t b + FTT01 12 22（4）設 n m 數(shù)據(jù)陣 E0 的秩為 r min(n -1, m) ，則存在r 個成分t1, t2 , tr ， 使得 E= t a+ t a +ETT01 1r rrF = t b + t b + FT+T 01 1r rr把tk = wk1x1 + wkm xm （ k = 1,2, r ），代入Y = t1b1 + tr br ，即得 p 個因變量的偏最小二乘回歸方程式 -533- y j = a j1x1 +

9、 a jm xm ，（ j = 1,2,m） （5）交叉有效性檢驗。 一般情況下，偏最小二乘法并不需要選用存在的 r 個成分t1 , t2 , tr 來建立回歸 式，而像主成分分析一樣，只選用前l(fā) 個成分（ l r ），即可得到預測能力較好的回歸模型。對于建模所需提取的主成分個數(shù)l ，可以通過交叉有效性檢驗來確定。 每次舍去第i 個觀測（ i = 1,2, n ），用余下的 n - 1 個觀測值按偏最小二乘回歸 方法建模，并考慮抽取h 個成分后擬合的回歸式，然后把舍去的第i 個觀測點代入所擬合的回歸方程式，得到 y j ( j = 1,2, p) 在第i 個觀測點上的預測值 y(i) j (h

10、) 。對 i = 1,2, n 重復以上的驗證，即得抽取 h 個成分時第 j 個因變量 y j ( j = 1,2, p)的 預測誤差平方和為 nPRESS j(h) = ( y i-j y (i) j (h)2i=1（ j = 1,2, p Y = ( y1 , y p)T 的預測誤差平方和為 pPRESS(h) = PRESS j (h) 。i=1另外，再采用所有的樣本點，擬合含h 個成分的回歸方程。這時，記第i 個樣本點的預測值為 yij (h) ，則可以定義 y j 的誤差平方和為 nSS j(h) = ( y i-j y i(jh)2i=1定義Y 的誤差平方和為pSS(h) = SS

11、 j (h)j=1當PRESS(h) 達到最小值時，對應的 h 即為所求的成分個數(shù)。通常，總有PRESS(h) 大于SS(h) ，而SS(h) 則小于SS(h -1) 。因此，在提取成分時，總希望比值PRESS(h) SS(h -1) 越小越好；一般可設定限制值為 0.05，即當-534- PRESS(h) SS(h -1) (1 - 0.05)2 = 0.952時，增加成分th 有利于模型精度的提高?；蛘叻催^來說，當 PRESS(h) SS(h -1) 0.952時，就認為增加新的成分th ，對減少方程的預測誤差無明顯的改善作用。 為此，定義交叉有效性為Q2 = 1 - PRESS(h) S

12、S(h -1)，這樣，在建模的每一步 h計算結束前，均進行交叉有效性檢驗，如果在第 h 步有Q2h1Q_h2(i)=1-press(i)/ss(i-1);elseQ_h2(1)=1;endif Q_h2(i)0.0975fprintf(提出的成分個數(shù) r=%d,i); r=i;breakendendbeta_z=t(:,1:r),ones(num,1)f0;%求Y 關于 t 的回歸系數(shù)beta_z(end,:)=;%刪除常數(shù)項xishu=w_star(:,1:r)*beta_z;每一列是一個回歸方程 %求Y 關于X 的回歸系數(shù)，且是針對標準數(shù)據(jù)的回歸系數(shù)， mu_x=mu(1:n);mu_y=

13、mu(n+1:end); sig_x=sig(1:n);sig_y=sig(n+1:end);-538- for i=1:mch0(i)=mu_y(i)-mu_x./sig_x*sig_y(i)*xishu(:,i);%計算原始數(shù)據(jù)的回歸方程的常數(shù) 項 endfor i=1:mxish(:,i)=xishu(:,i)./sig_x*sig_y(i);%計算原始數(shù)據(jù)的回歸方程的系數(shù)，每一列是一個回歸方程 endsol=ch0;xish常數(shù)項 %顯示回歸方程的系數(shù)，每一列是一個方程，每一列的第一個數(shù)是 save mydata x0 y0 num xishu ch0 xish計算得只要提出兩個成分t

14、1,2t 即可，交叉有效性Q2 = -0.1969 。 w 與 w* 的取值2hh見表 3，成分th 的得分th見表 4。 *表 3w 與 w 的取值hh表 4 成分th 的得分th標準化變量 y 關于成分t 的回歸模型如下k1y= r t + r t ， k = 1,2,3k1k 12k 2-539-No12345678910t1t2Not1t2-0.6429 -0.5914 110.3645 -0.7007 -0.7697 -0.1667 120.7433 -0.6983 -0.9074 0.5212 131.1867 0.757 0.6884 0.68 14-4.3898 0.76 -0

15、.4867 -1.1328 15-0.8232 -0.9738 -0.2291 0.0717 16-0.749 0.5211 -1.4037 0.0767 17-0.3929 0.2034 0.7436 0.2106 181.1993 -0.7827 1.7151 0.6549 191.0485 -0.3729 1.1626 -0.1668 201.9424 1.1294 自變量 w1w2*w1*w2x1-0.5899 -0.4688-0.5899 -0.3679x2-0.7713 0.5680-0.7713 0.6999x30.2389 0.67650.2389 0.6356 由于成分th

16、可以寫成原變量的標準化變量 x j 的函數(shù)，即有xxt = w x + w*+ w*h1h 12h 23h 3由此可得由成分t1 所建立的偏最小二乘回歸模型為yx= r (w x* + w x *+ w* x ) + r (w* x + w*+ w*x )k1k11 121231 32k12 122 232 3= (r+ r w* )x + (r+ r)x + (r+ r)xw*w*w*w*w*1k 112k 1211k 212k 2221k 312k323有關 rh = (rh1 , rh 2 , rh3 ) 的計算結果見表 5。表 5 回歸系數(shù) rh所以，有y= -0.0778x - 0.

17、4989x - 0.1322x1123yy= -0.1385x - 0.5244x - 0.0854x2123= -0.0604x - 0.1559x - 0.0073x3123將標準化變量 y , x (k = 1,2,3) 分別還原成原始變量 y , x (k = 1,2,3) ，則回歸方 kkkk程為 y1 = 47.0197 - 0.0167x1 - 0.8237x2 - 0.0969x3y2 = 612.5671 - 0.3509x1 - 10.2477x2 - 0.7412x3y3 = 183.9849 - 0.1253x1 - 2.4969x2 - 0.0518x3為了更直觀、迅

18、速地觀察各個自變量在解釋 yk (k = 1,2,3) 時的邊際作用，可以繪制回歸系數(shù)圖，見圖 1。這個圖是針對標準化數(shù)據(jù)的回歸方程的。 從回歸系數(shù)圖中可以立刻觀察到，腰圍變量在解釋三個回歸方程時起到了極為重要的作用。然而，與單杠及彎曲相比，跳高成績的回歸方程顯然不夠理想，三個自變量對它的解釋能力均很低。 -540- k123r10.34160.41610.1430r2-0.3364 -0.2908 -0.0652 為了考察這三個回歸方程的模型精度，我們以( yik , yik ) 為坐標值，對所有的樣本 點繪制預測圖。 yik 是第 k 個變量，第i 個樣本點( yik ) 的預測值。在這個預測圖上，如 果所有點都能在圖的對角線附近均勻分布，則方程的擬合值與原值差異很小，這個方程的擬合效果就是滿意的。體能訓練的預測圖見圖 2。 0-0.1-0.2-0.3-0.4-0.5-0.6單杠 彎曲 跳高 -0.7123圖 1回歸系數(shù)的直方圖2030015 200101005000510152001002003002502001501005000100200300圖 2 體能訓練預測圖 畫直方圖的 MATLAB 程序為：bar(xishu)畫體能訓練的預測圖的 MATLAB 程序如下：load mydata numch0=repmat(ch0,num,1); yhat=ch0+x0*xis