偏最小二乘法

數學優化技術之一

偏最小二乘法是一種數學優化技術,它通過最小化誤差的平方和找到一組數據的最佳函數匹配。用最簡的方法求得一些絕對不可知的真值,而令誤差平方之和為最小。很多其他的優化問題也可通過最小化能量或最大化熵用最小二乘形式表達。

基本公式


偏最小二乘回歸≈多元線性回歸分析+典型相關分析+主成分分析

基本特點


與傳統多元線性回歸模型相比,偏最小二乘回歸的特點是:
(1)能夠在自變數存在嚴重多重相關性的條件下進行回歸建模;
(2)允許在樣本點個數少於變數個數的條件下進行回歸建模;
(3)偏最小二乘回歸在最終模型中將包含原有的所有自變數;
偏最小二乘法
偏最小二乘法
(4)偏最小二乘回歸模型更易於辨識系統信息與雜訊(甚至一些非隨機性的雜訊);
(5)在偏最小二乘回歸模型中,每一個自變數的回歸係數將更容易解釋。
在計算方差和協方差時,求和號前面的係數有兩種取法:當樣本點集合是隨機抽取得到時,應該取1/(n-1);如果不是隨機抽取的,這個係數可取1/n。

診斷方法


經驗式診斷法
1、在自變數的簡單相關係數矩陣中,有某些自變數的相關係數值較大。
2、回歸係數的代數符號與專業知識或一般經驗相反;或者,它同該自變數與y的簡單相關係數符號相反。
3、對重要自變數的回歸係數進行t檢驗,其結果不顯著。
特別典型的是,當F檢驗能在高精度下通過,測定係數R2的值亦很大,但自變數的t檢驗卻全都不顯著,這時,多重相關性的可能性將很大。
4、如果增加(或刪除)一個變數,或者增加(或刪除)一個觀測值,回歸係數的估計值發生了很大的變化。
5、重要自變數的回歸係數置信區間明顯過大。
6、在自變數中,某一個自變數是另一部分自變數的完全或近似完全的線性組合。
7、對於一般的觀測數據,如果樣本點的個數過少,樣本數據中的多重相關性是經常存在的。
但是,採用經驗式方法診斷自變數系統中是否確實存在多重相關性,並不十分可靠,另一種較正規的方法是利用統計檢驗(回歸分析),檢查每一個自變數相對其它自變數是否存在線性關係
方差膨脹因子診斷法
最常用的多重相關性的正規診斷方法是使用方差膨脹因子。自變數xj的方差膨脹因子記為(VIF)j,它的計算方法為
(4-5)(VIF)j=(1-Rj2)-1
式中,Rj2是以xj為因變數時對其它自變數回歸的複測定係數。
所有xj變數中最大的(VIF)j通常被用來作為測量多重相關性的指標。一般認為,如果最大的(VIF)j超過10,常常表示多重相關性將嚴重影響最小二乘的估計值。
(VIF)j被稱為方差膨脹因子的原因,是由於它還可以度量回歸係數的估計方差與自變數線性無關時相比,增加了多少。
不妨假設x1,x2,…,xp均是標準化變數。採用最小二乘法得到回歸係數向量B,它的精度是用它的方差來測量的。B的協方差矩陣為
Cov(B)=σ2(X'X)-1
式中,σ2是誤差項方差。所以,對於回歸係數bj,有
Var(bj)=σ2cjj
cjj是(X'X)-1矩陣中第j個對角元素。可以證明,
cjj=(VIF)j

嶺回歸分析


嶺回歸估計量簡介
嶺回歸分析是一種修正的最小二乘估計法,當自變數系統中存在多重相關性時,它可以提供一個比最小二乘法更為穩定的估計,並且回歸係數的標準差也比最小二乘估計的要小。
根據高斯——馬爾科夫定理,多重相關性並不影響最小二乘估計量的無偏性和最小方差性。但是,雖然最小二乘估計量在所有線性無偏估計量中是方差最小的,但是這個方差卻不一定小。於是可以找一個有偏估計量,這個估計量雖然有微小的偏差,但它的精度卻能夠大大高於無偏的估計量。
在應用嶺回歸分析時,它的計算大多從標準化數據出發。對於標準化變數,最小二乘的正規方程為
rXXb=ryX
式中,rXX是X的相關係數矩陣,ryX是y與所有自變數的相關係數向量。
嶺回歸估計量是通過在正規方程中引入有偏常數c(c≥0)而求得的。它的正規方程為+
(4-8)(rXX+cI)bR=ryX
所以,在嶺回歸分析中,標準化回歸係數為
(4-9)bR=(rXX+cI)-1ryX
嶺回歸估計量性質
(1)嶺回歸係數是一般最小二乘準則下回歸係數的線性組合,即
(4-10)bR=(I+crXX-1)-1b
(2)記β是總體參數的理論值。當β≠0時,可以證明一定存在一個正數c0,使得當0
(4-11)E||bR-β||2≤E||b-β||2
(3)嶺回歸估計量的絕對值常比普通最小二乘估計量的絕對值小,即
(4-12)||bR||<||b||
嶺回歸估計量的質量取決於偏倚係數c的選取。c的選取不宜過大,因為
E(bR)=(I+crXX-1)-1E(b)=(I+crXX-1)-1β
關於偏倚係數c的選取尚沒有正規的決策準則,目前主要以嶺跡和方差膨脹因子為依據。嶺跡是指p-1個嶺回歸係數估計量對不同的c值所描繪的曲線(c值一般在0~1之間)。在通過檢查嶺跡和方差膨脹因子來選擇c值時,其判斷方法是選擇一個儘可能小的c值,在這個較小的c值上,嶺跡中的回歸係數已變得比較穩定,並且方差膨脹因子也變得足夠小。
從理論上,最佳的c值是存在的,它可以使估計量的偏差和方差的組合效應達到一個最佳水準。然而,困難卻在於c的最優值對不同的應用而有所不同,對其選擇還只能憑經驗判斷。

其他補救方法


去掉相關性變數
由於變數間多重相關性的形式十分複雜,而且還缺乏十分可靠的檢驗方法,刪除部分多重相關變數的做法常導致增大模型的解釋誤差,將本應保留的系統信息捨棄,使得接受一個錯誤結論的可能和做出錯誤決策的風險都不斷增長。另一方面,在一些經濟模型中,從經濟理論上要求一些重要的解釋變數必須被包括在模型中,而這些變數又存在多重相關性。這時採用剔除部分相關變數的做法就不符合實際工作的要求。
然而,在實際工作中,由於時間、經費以及客觀條件的限制,增大樣本容量的方法常常是不可行的。
變數轉換
此法可用來削弱多重相關性的嚴重性。一階差分回歸模型有可能減少多重相關性的嚴重性。然而,一階差分變換又帶來了一些其它問題。差分后的誤差項可能不滿足總體模型中關於誤差項不是序列相關的假定。事實上,在大部分情形下,在原來的誤差項是不自相關的條件下,一階差分所得到的誤差項將會是序列相關的。而且,由於差分方法損失了一個觀察值,這在小樣本的情況下是極不可取的。另外,一階差分方法在截面樣本中是不宜利用的。
1主成分分析
主成分分析的計算結果必然受到重疊信息的影響。因此,當人為地採用一些無益的相關變數時,無論從方向上還是從數量上,都會扭曲客觀結論。在主成分分析之前,對變數系統的確定必須是慎之又慎的。
2特異點的發現
第i個樣本點(樣本量為n)對第h主成分的貢獻率是
(5-32)CTR(i)=Fh2(i)/(nλh)(若遠超過1/n,為特異點)
3典型相關分析
從某種意義上說,多元回歸分析、判別分析或對應分析等許多重要的數據分析方法,都可以歸結為典型相關分析的一種特例,同時它還是偏最小二乘回歸分析的理論基石。
典型相關分析,是從變數組X中提取一個典型成分F=Xa,再從變數組Y中提取一個成分G=Yb,在提取過程中,要求F與G的相關程度達到最大。
在典型相關分析中,採用下述原則尋優,即
max=aX'Yba'X'Xa=1,b'Y'Yb=1
其結果為,a是對應於矩陣V11-1V12V22-1V21最大特徵值的特徵向量,而b是對應於矩陣V22-1V21V11-1V12最大特徵值的特徵向量,這兩個最大特徵值相同。其中,
V11=X'X,V12=X'Y,V22=Y'Y。
F與G之間存在著明顯的換算關係。
有時只有一個典型成分還不夠,還可以考慮第二個典型成分。
多因變數的偏最小二乘回歸模型

建模方法


設有q個因變數和p個自變數。為了研究因變數與自變數的統計關係,觀測了n個樣本點,由此構成了自變數與因變數的數據表X和Y。偏最小二乘回歸分別在X與Y中提取出t和u,要求:(1)t和u應儘可能大地攜帶它們各自數據表中的變異信息;(2)t和u的相關程度能夠達到最大。在第一個成分被提取后,偏最小二乘回歸分別實施X對t的回歸以及Y對t的回歸。如果回歸方程已經達到滿意的精度,則演演算法終止;否則,將利用X被t解釋后的殘餘信息以及Y被t解釋后的殘餘信息進行第二輪的成分提取。如此往複,直到能達到一個較滿意的精度為止。若最終對X共提取了多個成分,偏最小二乘回歸將通過施行yk對X的這些成分的回歸,然後再表達成yk關於原自變數的回歸方程。

計算方法


計算方法1

首先將數據做標準化處理。X經標準化處理后的數據矩陣記為E0=(E01,…,E0p)n×p,Y的相應矩陣記為F0=(F01,…,F0q)n×q。
第一步記t1是E0的第一個成分,t1=E0w1,w1是E0的第一個軸,它是一個單位向量,即||w1||=1。
記u1是F0的第一個成分,u1=F0c1,c1是F0的第一個軸,並且||c1||=1。
於是,要求解下列優化問題,即
(7-1)
記θ1=w1'E0'F0c1,即正是優化問題的目標函數值。
採用拉格朗日演演算法,可得
(7-8)E0'F0F0'E0w1=θ12w1
(7-9)F0'E0E0'F0c1=θ12c1
所以,w1是對應於E0'F0F0'E0矩陣最大特徵值的單位特徵向量,而c1是對應於F0'E0E0'F0矩陣最大特徵值θ12的單位特徵向量。
求得軸w1和c1后,即可得到成分
t1=E0w1
u1=F0c1
然後,分別求E0和F0對t1的回歸方程
(7-10)E0=t1p1'+E1
(7-12)F0=t1r1'+F1
式中,回歸係數向量是
(7-13)p1=E0't1/||t1||2
(7-15)r1=F0't1/||t1||2
而E1和F1分別是兩個方程的殘差矩陣。
第二步用殘差矩陣E1和F1取代E0和F0,然後,求第二個軸w2和c2以及第二個成分t2,u2,有
t2=E1w2
u2=F1c2
θ2==w2'E1'F1c2
w2是對應於E1'F1F1'E1矩陣最大特徵值的單位特徵向量,而c2是對應於F1'E1E1'F1矩陣最大特徵值θ22的單位特徵向量。計算回歸係數
p2=E1't2/||t2||2
r2=F1't2/||t2||2
因此,有回歸方程
E1=t2p2'+E2
F1=t2r2'+F2
如此計算下去,如果X的秩是A,則會有
(7-16)E0=t1p1'+…+tApA'
(7-17)F0=t1r1'+…+tArA'+FA
由於t1,…,tA均可以表示成E01,…,E0p的線性組合,因此,式(7-17)還可以還原成yk*=F0k關於xj*=E0j的回歸方程形式,即
yk*=αk1x1*+…+αkpxp*+FAk,k=1,2,…,q
FAk是殘差矩陣FA的第k列。
3交叉有效性
如果多一個成分而少一個樣本的預測誤差平方和(所有因變數和預測樣本相加)除以少一個成分的誤差平方和(所有的因變數和樣本相加)小於0.952,則多一個成分是值得的。

計算方法2

用下述原則提取自變數中的成分t1,是與原則式(7-1)的結果完全等價的,即
(7-24)
(1)求矩陣E0'F0F0'E0最大特徵值所對應的單位特徵向量w1,求成分t1,得
t1=E0w1
E1=E0-t1p1'
式中,p1=E0't1/||t1||2
(2)求矩陣E1'F0F0'E1最大特徵值所對應的單位特徵向量w2,求成分t2,得
t2=E1w2
E2=E1-t2p2'
式中,p2=E1't2/||t2||2
……
(m)至第m步,求成分tm=Em-1wm,wm是矩陣Em-1'F0F0'Em-1最大特徵值所對應的單位特徵向量.
如果根據交叉有效性,確定共抽取m個成分t1,…,tm可以得到一個滿意的觀測模型,則求F0在t1,…,tm上的普通最小二乘回歸方程為
F0=t1r1'+…+tmrm'+Fm
偏最小二乘回歸的輔助分析技術
1精度分析
定義自變數成分th的各種解釋能力如下
(1)th對某自變數xj的解釋能力
(8-1)Rd(xj;th)=r2(xj,th)
(2)th對X的解釋能力
(8-2)Rd(X;th)=[r2(x1,th)+…+r2(xp,th)]/p
(3)t1,…,tm對X的累計解釋能力
(8-3)Rd(X;t1,…,tm)=Rd(X;t1)+…+Rd(X;tm)
(4)t1,…,tm對某自變數xj的累計解釋能力
(8-4)Rd(xj;t1,…,tm)=Rd(xj;t1)+…+Rd(xj;tm)
(5)th對某因變數yk的解釋能力
(8-5)Rd(yk;th)=r2(yk,th)
(6)th對Y的解釋能力
(8-6)Rd(Y;th)=[r2(y1,th)+…+r2(yq,th)]/q
(7)t1,…,tm對Y的累計解釋能力
(8-7)Rd(Y;t1,…,tm)=Rd(Y;t1)+…+Rd(Y;tm)
(8)t1,…,tm對某因變數yk的累計解釋能力
(8-8)Rd(yk;t1,…,tm)=Rd(yk;t1)+…+Rd(yk;tm)
2自變數xj在解釋因變數集合Y的作用
xj在解釋Y時作用的重要性,可以用變數投影重要性指標VIPj來測度
VIPj2=p[Rd(Y;t1)w1j2+…+Rd(Y;tm)wmj2]/[Rd(Y;t1)+…+Rd(Y;tm)]
式中,whj是軸wh的第j個分量。注意VIP12+…+VIPp2=p
3特異點的發現
定義第i個樣本點對第h成分th的貢獻率Thi2,用它來發現樣本點集合中的特異點,即
(8-10)Thi2=thi2/((n-1)sh2)
式中,sh2是成分th的方差。
由此,還可以測算樣本點i對成分t1,…,tm的累計貢獻率
(8-11)Ti2=T1i2+…+Tmi2
Ti2≥m(n2-1)F0.05(m,n-m)/(n2(n-m))
時,可以認為在95%的檢驗水平上,樣本點i對成分t1,…,tm的貢獻過大。
單因變數的偏最小二乘回歸模型
1簡化演演算法
第一步已知數據E0,F0,由於u1=F0,可得
w1=E0'F0/||E0'F0||
t1=E0w1
p1=E0't1/||t1||2
E1=E0-t1p1'
檢驗交叉有效性。若有效,繼續計算;否則只提取一個成分t1。
第h步(h=2,…,m)已知數據Eh-1,F0,有
wh=Eh-1'F0/||Eh-1'F0||
th=Eh-1wh
ph=Eh-1'th/||th||2
Eh=Eh-1-thph'
檢驗交叉有效性。若有效,繼續計算h+1步;否則停止求成分的計算。
這時,得到m個成分t1,…,tm,實施F0在t1,…,tm上的回歸,得
F0^=r1t1+…+rmtm
由於t1,…,tm均是E0的線性組合,即
th=Eh-1wh=E0wh*
所以F0^可寫成E0的線性組合形式,即
F0^=r1E0w1*+…+rmE0wm*=E0[r1w1*+…+rmwm*]
最後,也可以變換成y對x1,…,xp的回歸方程
y^=α0+α1x1+…+αpxp
  • 目錄