【建模應用】PLS偏最小二乘回歸原理與應用

本文轉載自查看原文 2018-05-13 11:13 42066 matlab/ 算法/ 數學建模

出處：http://www.cnblogs.com/duye/p/9031511.html

1.回歸

　　”回歸“一詞來源於對父母身高對於子女身高影響的研究。有人對父母的身高與子女身高做統計，發現除了父母高則子女普遍高的常識性結論外，子女的身高總是“趨向”於人類平均身高，最早“回歸”一詞就來源於此，即子女的身高總是回歸於人類平均身高。
　　現代意義上的回歸，是研究因變量對自變量的依賴關系的一種統計分析方法，目的是通過自變量的給定值來估計或預測因變量的均值。它可用於預測、時間序列建模以及發現各種變量之間的因果關系。簡單地說，回歸就是去分析因變量與自變量之間的關系，從而為分析數據、預測數據提供科學的、合理的方法。

2.回歸的方法

目前常用的回歸方法有以下幾種：

多元線性回歸：入門練習常見，但用在模型上基本無用，因為多個變量之間難免有復雜的相關性，多元線性回歸不能處理多個自變量之間的“糅合”性。
邏輯回歸：當預測的是1/0時使用。這也是機器學習中的一種常用二分類方法。
主成分回歸：類似於主成分分析，將具有相關性的多維變量降維到互不相關的少數幾維變量上，可以處理變量之間具有相關性的情況。

嶺回歸：同上，但方法不同。

偏最小二乘回歸：當數據量小，甚至比變量維數還小，而相關性又比較大時使用，這個方法甚至優於主成分回歸。

3.回歸的檢驗

　　建模是最重要的，但好的回歸模型是需要檢驗的，否則你的模型就會顯得“蒼白無力”。回歸結果好與壞，應該怎么檢驗呢？從以下幾方入手:

自變量與因變量是否具有預期的關系。如果有非常不符合邏輯的系數，我們就應該考慮剔除它了。
自變量對模型是否有幫助。如果自變量的系數為零（或非常接近零），我們認為這個自變量對模型沒有幫助，統計檢驗就用來計算系數為零的概率。如果統計檢驗返回一個小概率值（p值），則表示系數為零的概率很小。如果概率小於0.05，匯總報告上概率（Probability）旁邊的一個星號（*）表示相關自變量對模型非常重要。換句話說，其系數在95%置信度上具有統計顯著性。
殘差是否有空間聚類。殘差在空間上應該是隨機分布的，而不應該出現聚類。
模型是否出現了傾向性。如果我們正確的構建了回歸分析模型，那么模型的殘差會符合完美的正態分布，其圖形為鍾形曲線。
自變量中是否存在冗余。建模的過程中，應盡量去選擇表示各個不同方面的自變量，也就是盡量避免傳達相同或相似信息的自變量。 評估模型的性能。評估矯R2值，有時還要加上修正的Akaike信息准則/Akaike's information criterion (AIC)，效果是否好。

偏最小二乘回歸

1.偏最小二乘回歸的思想：

　　一般來說，能用主成分分析就能用偏最小二乘。偏最小二乘集成了主成分分析、典型相關分析、線性回歸分析的優點。在普通多元線形回歸的應用中，我們常受到許多限制。最典型的問題就是：自變量之間的多重相關性。並且有的時候樣例很少，甚至比變量的維度還少，變量之間又存在多重相關性。偏最小二乘回歸就是為解決這些棘手的問題而生的。
　　舉個例子，比如現在，有一堆因素（X1,X2,...Xn）(這是自變量)，這些因素可以導致（Y1,Y2,...Yn）（這是因變量），給的樣例很少，而我們又完全不清楚自變量之間、因變量之間存在的關系，這時問自變量與因變量之間到底是一個什么關系？這就是偏最小二乘要解決的問題。

2.偏最小二乘回歸建模原理：

　　設有 q個因變量{y1,…,yq}和p自變量{x1,…,xp}。為了研究因變量和自變量的統計關系,我們觀測了n個樣本點,由此構成了自變量與因變量的數據表X={x1,…,xp}和.Y={y1,…,yp}。偏最小二乘回歸分別在X與Y中提取出成分t1和u1(也就是說, t1是x1,x2,…,xq的線形組合, u1是y1,y2,…,yp的線形組合).在提取這兩個成分時,為了回歸分析的需要,有下列兩個要求:

　　(1) t1和u1應盡可能大地攜帶他們各自數據表中的變異信息;

　　(2) t1與u1的相關程度能夠達到最大。

　　這兩個要求表明，t1和u1應盡可能好的代表數據表X和Y,同時自變量的成分t1對因變量的成分u1又有最強的解釋能力。

　　在第一個成分t1和u1被提取后，偏最小二乘回歸分別實施X對t1的回歸以及 Y對u1的回歸。如果回歸方程已經達到滿意的精度，則算法終止；否則,將利用 X被t1解釋后的殘余信息以及Y 被t2解釋后的殘余信息進行第二輪的成分提取。如此往復，直到能達到一個較滿意的精度為止。若最終對 X共提取了 m個成分t1,t2,…,tm，偏最小二乘回歸將通過實施yk對t1,t2,…,tm的回歸,然后再表達成yk關於原變量X1,X2,…,Xq的回歸方程,k=1,2,…,p。

3.推導偏最小二乘回歸：

　　為了徹底理解偏最小二乘回歸，我建議下面的步驟你都親自推導一遍。相信經過下面的推導，能讓你對偏最小二乘有一個更加清晰的認識。

step1:數據說明與標准化

　　數據矩陣E0，F0，其中E0為自變量矩陣，每一行是一個樣例，每一列代表了一個維度的變量；F0是因變量矩陣，解釋同E0。

　　數據標准化即，要將數據中心化，方法是每個樣本都做如下操作：減去一個維度變量的均值除以該維度的標准差。以下設E0，F0都為標准化了的數據。即：自變量經標准化處理后的數據矩陣記為E0（nm），因變量經標准化處理后的數據矩陣記為F0（np）。

step2:求符合要求的主成分（☆）

　　即求自變量與因變量的第一對主成分t1和u1，根據主成分原理，要求t1與u1的方差達到最大，這是因為：方差最大則表示的信息就越多。另一方面，又要求t1對u1有最大的解釋能力，由典型相關分析的思路知，t1與u1的相關度達到最大值。

　　因此，綜合上述兩點，我們只要要求t1與u1的協方差達到最大，即：

　　　　　　　　　　Cov(t1,u1) –> max

　　而且，t1是X的線性組合，那么權重系數設為W1，即t1 = E0W1，同理，u1是Y的線性組合，u1 = F0C1。同時又要求，W1與C1同為單位向量，問題的數學表達式為：

　　　　　　　　　　max <E0w1,F0c1>

　　　　　　　　　　S.T.

　　　　　　　　　　　　||W1|| = 1;

　　　　　　　　　　　　||c1|| = 1

　　這就是一個條件極值的問題，你可以采用拉格朗日方法求解（如果你還有興趣，可以查閱高數課本，當然，你也可以直接看結論，這里我只給出結論，推導省略，實際上推導並不影響你理解）。

　　通過拉格朗日求解，知w1就是矩陣E0‘F0F0’E0的對應於最大特征值的特征向量，c1就是矩陣F0’E0E0’F0對應於最大特征值的最大特征向量，均單位化。

　　有了權系數w1，c1，自然可以求得主成分t1，u1。至此，第一對主成分完成。

step3:建立主成分與原自變量、因變量之間的回歸（☆）

　　建立E0，F0對t1，u1的三個回歸方程，如下：

　　式中，回歸系數向量是：

　　而E1,F*1,F1分別是三個回歸方程的殘差矩陣.

step4:繼續求主成分，直到滿足要求

　　用殘差矩陣E1和F1取代E0和F0,然后,求第二個軸w2和c2以及第二個成分t2,u2,有

　　重新執行step3。直到求出所有主成分或者滿足要求（后面說明）。

step5:推導因變量之於自變量的回歸表達式（☆）

　　如此經過step3-step4反復，若E0的秩為A，則可以求出：

　　由於t1….tA都可以表示E01，E02….E0q的線性組合，那么就自然還原成下面的形式：

　　Fak為殘差矩陣Fa的第k列。這樣，就求出了回歸方程。

step6:檢驗-交叉有效性（☆）

　　這是最后一步，也是非常重要的一步。下面要討論的問題是在現有的數據表下,如何確定更好的回歸方程。在許多情形下,偏最小二乘回歸方程並不需要選用全部的成分進行回歸建模,而是可以象在主成分分析一樣,采用截尾的方式選擇前m 個成分,僅用這m 個后續的成分就可以得到一個預測性較好的模型。事實上,如果后續的成分已經不能為解釋因變量提供更有意義的信息時,采用過多的成分只會破壞對統計趨勢的認識,引導錯誤的預測結論。

　　下面的問題是怎樣來確定所應提取的成分個數。

　　在偏最小二乘回歸建模中,究竟應該選取多少個成分為宜,這可通過考察增加一個新的成分后,能否對模型的預測功能有明顯的改進來考慮。采用類似於抽樣測試法的工作方式,把所有n個樣本點分成兩部分:第一部分除去某個樣本點i的所有樣本點集合(共含n-1個樣本點),用這部分樣本點並使用h個成分擬合一個回歸方程;第二部分是把剛才被排除的樣本點i代入前面擬合的回歸方程,得到yj在樣本點i上的擬合值。對於每一個i=1,2,…,n,重復上述測試,則可以定義yj的預測誤差平方和為PRESShj。有：

　　定義Y 的預測誤差平方和為PRESSh,有

　　顯然,如果回歸方程的穩健性不好,誤差就很大,它對樣本點的變動就會十分敏感,這種擾動誤差的作用,就會加大PRESSh的值。

　　另外,再采用所有的樣本點,擬合含h 個成分的回歸方程。這是,記第i個樣本點的預測值為,則可以記yj的誤差平方和為SShj,有

　　定義Y的誤差平方和為SSh,有

　　定義稱為交叉有效性，對於每一個變量yk，定義

　　對於全部因變量Y,成分th交叉有效性定義為

　　用交叉有效性測量成分th對預測模型精度的邊際貢獻有如下兩個尺度。

　　(1) 當時, th成分的邊際貢獻是顯著的。顯而易見,是完全等價的決策原則。

　　(2) 對於k=1,2,…,q,至少有一個k,使得。這時增加成分th,至少使一個因變量yk的預測模型得到顯著的改善,因此,也可以考慮增加成分th是明顯有益的。

實現偏最小二乘回歸算法步驟：

　　上面推導了偏最小二乘回歸，分析了其中的原理。為了使得在實際應用中更加快速的使用偏最小二乘回歸，在此，貼上實現偏最小二乘法實現的簡潔步驟，需說明的是，下面算法來自司守奎老師《數學建模算法與應用》一書，該書推導過程跨度大，個人認為不適合新手直接閱讀，建議你在理解了上述第二部分后再去閱讀此書“偏最小二乘回歸”章節，定會有更加高層次的認識。步驟如下：

四、MATLAB實例以及實現

　　有必要貼出偏最小二乘的簡單建模應用，並用matlab去是實現之，你可以按照上述步驟，通過基本的運算如求矩陣特征值等，來實現，也可以使用matlab工具箱方法實現之，下面給出的依舊是一個來自司守奎老師書本上的案例：

　　例：采用蘭納胡德（Linnerud）給出的關於體能訓練的數據進行偏小二乘回歸建模。在這個數據系統中被測的樣本點，是某健身俱樂部的 20 位中年男子。被測變量分為兩組。第一組是身體特征指標 X ，包括：體重、腰圍、脈搏。第二組變量是訓練結果指標Y ，包括：單杠、彎曲、跳高。原始數據見表 1。表 2 給出了這 6 個變量的簡單相關系數矩陣。從相關系數矩陣可以看出，體重與腰圍是正相關的；體重、腰圍與脈搏負相關；而在單杠、彎曲與跳高之間是正相關的。從兩組變量間的關系看，單杠、彎曲和跳高的訓練成績與體重、腰圍負相關，與脈搏正相關。

　　可以利用如下的MATLAB程序：

 1 clc,clear  2 load pz.txt %原始數據存放在純文本文件 pz.txt 中  3 mu=mean(pz);sig=std(pz); %求均值和標准差  4 rr=corrcoef(pz); %求相關系數矩陣  5 data=zscore(pz); %數據標准化,變量記做 X*和 Y*
 6 n=3;m=3; %n 是自變量的個數,m 是因變量的個數  7 x0=pz(:,1:n);y0=pz(:,n+1:end); %原始的自變量和因變量數據  8 e0=data(:,1:n);f0=data(:,n+1:end); %標准化后的自變量和因變量數據  9 -679-
10 num=size(e0,1);%求樣本點的個數 11 chg=eye(n); %w 到 w*變換矩陣的初始化 12 for i=1:n 13 %以下計算 w，w*和 t 的得分向量， 14 matrix=e0'*f0*f0'*e0; 15 [vec,val]=eig(matrix); %求特征值和特征向量 16 val=diag(val); %提出對角線元素，即提出特征值 17 [val,ind]=sort(val,'descend'); 18 w(:,i)=vec(:,ind(1)); %提出最大特征值對應的特征向量 19 w_star(:,i)=chg*w(:,i); %計算 w*的取值 20 t(:,i)=e0*w(:,i); %計算成分 ti 的得分 21 alpha=e0'*t(:,i)/(t(:,i)'*t(:,i)); %計算 alpha_i 22 chg=chg*(eye(n)-w(:,i)*alpha'); %計算 w 到 w*的變換矩陣
23 e=e0-t(:,i)*alpha'; %計算殘差矩陣
24 e0=e; 25 %以下計算 ss(i)的值 26 beta=t\f0; %求回歸方程的系數，數據標准化，沒有常數項 27 cancha=f0-t*beta; %求殘差矩陣 28 ss(i)=sum(sum(cancha.^2)); %求誤差平方和 29 %以下計算 press(i) 30 for j=1:num 31 t1=t(:,1:i);f1=f0; 32 she_t=t1(j,:);she_f=f1(j,:); %把舍去的第 j 個樣本點保存起來 33 t1(j,:)=[];f1(j,:)=[]; %刪除第 j 個觀測值 34 beta1=[t1,ones(num-1,1)]\f1; %求回歸分析的系數,這里帶有常數項 35 cancha=she_f-she_t*beta1(1:end-1,:)-beta1(end,:); %求殘差向量 36 press_i(j)=sum(cancha.^2); %求誤差平方和 37 end 38 press(i)=sum(press_i); 39 Q_h2(1)=1; 40 if i>1, Q_h2(i)=1-press(i)/ss(i-1); end 41 if Q_h2(i)<0.0975
42 fprintf('提出的成分個數 r=%d',i); break
43 end 44 end 45 beta_z=t\f0; %求 Y*關於 t 的回歸系數 46 xishu=w_star*beta_z; %求 Y*關於 X*的回歸系數，每一列是一個回歸方程 47 mu_x=mu(1:n);mu_y=mu(n+1:end); %提出自變量和因變量的均值 48 sig_x=sig(1:n);sig_y=sig(n+1:end); %提出自變量和因變量的標准差 49 ch0=mu_y-(mu_x./sig_x*xishu).*sig_y; %計算原始數據回歸方程的常數項 50 for i=1:m 51 xish(:,i)=xishu(:,i)./sig_x'*sig_y(i); %計算原始數據回歸方程的系數
52 end 53 sol=[ch0;xish] %顯示回歸方程的系數，每一列是一個方程，每一列的第一個數是常數項 54 save mydata x0 y0 num xishu ch0 xish