梁飛豹,林同華

(福州大學(xué) 數(shù)學(xué)與計(jì)算機(jī)科學(xué)學(xué)院,福建 福州　350108)

?

線性模型有偏估計(jì)的一種新算法

梁飛豹,林同華

(福州大學(xué) 數(shù)學(xué)與計(jì)算機(jī)科學(xué)學(xué)院,福建 福州350108)

摘要:文章以降低預(yù)測(cè)殘差平方和為目標(biāo),基于嶺估計(jì)增大回歸系數(shù)矩陣的對(duì)角元素的思想,提出一種利用高斯消去變換工具的線性模型參數(shù)估計(jì)法,并進(jìn)行數(shù)據(jù)模擬實(shí)驗(yàn),最后通過(guò)平均預(yù)測(cè)殘差平方和以及平均殘差平方和的箱線圖來(lái)對(duì)比新算法和最小二乘估計(jì)及嶺估計(jì)的優(yōu)良性,說(shuō)明滿足一定條件時(shí),新算法在估計(jì)精度和穩(wěn)定性上優(yōu)于這2種方法。

關(guān)鍵詞:預(yù)測(cè)殘差平方和;參數(shù)估計(jì);嶺估計(jì);高斯消去;箱線圖

考慮線性模型(Y,Xβ,σ2I):

(1)

為了減小線性模型的奇異性問(wèn)題,研究者引入有偏估計(jì)，常見(jiàn)的有偏估計(jì)有壓縮估計(jì)、Liu估計(jì)、主成分估計(jì)、嶺估計(jì)等。其中,壓縮估計(jì)[3]是對(duì)LS估計(jì)的均勻壓縮,它解決了因設(shè)計(jì)矩陣復(fù)共線性強(qiáng)而得到絕對(duì)值太大的參數(shù)估計(jì)的問(wèn)題;但所含的表達(dá)式包含未知參數(shù)β,實(shí)際過(guò)程中對(duì)它的估計(jì)存在一定的難度,并且在以MSE(mean squared error)作為估計(jì)的優(yōu)良性準(zhǔn)則時(shí),壓縮估計(jì)優(yōu)于LS需要回歸因子的復(fù)共線性不能太強(qiáng)。Liu估計(jì)[4]是關(guān)于LS估計(jì)的線性變換,但它的估計(jì)取值依賴于最小二乘估計(jì)。主成分估計(jì)[5]在減少對(duì)系統(tǒng)影響較小的自變量個(gè)數(shù)的同時(shí)也降低了與參數(shù)真值間的均方誤差,是一種線性約束下的最小二乘的解;但在主成分估計(jì)的過(guò)程中,XTX的特征值λi?0說(shuō)法較籠統(tǒng),如何保留主成分的個(gè)數(shù)是個(gè)重要的問(wèn)題。嶺估計(jì)[6]是在設(shè)計(jì)陣計(jì)算中引入偏參數(shù)k,通過(guò)對(duì)嶺參數(shù)k的合理取值有效地改善復(fù)共線性所帶來(lái)的病態(tài)性;而偏參數(shù)k擴(kuò)展成偏參數(shù)矩陣K,又得到廣義嶺估計(jì)。關(guān)于嶺參數(shù)的選擇,方法眾多[6-8]。在實(shí)際工作中,嶺估計(jì)是運(yùn)用較為廣泛的一種有偏估計(jì),許多研究者對(duì)嶺估計(jì)做了不同程度的改進(jìn),以期望縮小均方誤差、提高精度。2007年,Kaciranlar將壓縮估計(jì)和Swindel提出的修正嶺估計(jì)相結(jié)合,提出了兩參數(shù)估計(jì)[9],它是最小二乘估計(jì)、嶺估計(jì)、Liu估計(jì)、壓縮估計(jì)的綜合。2010年,楊虎和常新峰綜合嶺估計(jì)和Liu估計(jì)提出另一種兩參數(shù)估計(jì),并在MSE準(zhǔn)則下證明了該估計(jì)的優(yōu)良性[10]。

現(xiàn)有的有偏估計(jì)還難以找到一個(gè)最好的估計(jì),使得這個(gè)估計(jì)在有偏估計(jì)類中關(guān)于均方誤差最小。本文從近似奇異矩陣的求逆問(wèn)題角度出發(fā),提出一種新的線性模型有偏估計(jì)。

1高斯消去變換

(2)

性質(zhì)1Tkk(Tkk(A))=A,即對(duì)A連續(xù)施行2次消去變換,其結(jié)果都是A不變。

性質(zhì)2Tkk(Tll(A))=Tll(Tkk(A)),即對(duì)A施行2次消去變換,順序交換,其結(jié)果相同。

性質(zhì)3Tkk{Tll[Tmm(A)]}=Tmm{Tll[Tkk(A)]},即對(duì)A施行3次消去變換,結(jié)果與順序無(wú)關(guān)。

性質(zhì)4A-1=Tpp{T(p-1)(p-1){…T11(A)}},即按順序?qū)個(gè)對(duì)角元素作消去變換,得到逆矩陣。

對(duì)于線性模型(1)式,觀察增廣矩陣A=(XTXXTY),當(dāng)XTX可逆時(shí),對(duì)增廣矩陣A施行消去變換T11,T22,…,Tpp后,由(2)式可得:

可以看到此時(shí)矩陣B的最后一列即為線性模型(1)式的最小二乘估計(jì)。

由(2)式可知,在每次消去變換時(shí)主元充當(dāng)分母,如果主元太接近于0,計(jì)算機(jī)運(yùn)算過(guò)程中將不可避免地產(chǎn)生舍入誤差,其累積將產(chǎn)生較大的相對(duì)誤差。因此對(duì)增廣矩陣A施行消去變換,應(yīng)盡量避免選擇太小的主元。例如,每次選擇絕對(duì)值最大的未作變換的對(duì)角元素作為主元。

(3)

(4)

因此,LS估計(jì)的消去變換其實(shí)就是逐漸添加變量因子與觀測(cè)向量Y構(gòu)成線性方程并求解回歸系數(shù)的過(guò)程,沒(méi)有考慮逐漸添加的變量因子與Y的線性關(guān)系而僅以對(duì)角元素的大小作為主元選取的原則,有失客觀?？紤]到自變量因子與觀測(cè)向量的線性關(guān)系,本文采用逐步添加對(duì)觀測(cè)向量方差貢獻(xiàn)率最大的自變量因子與之構(gòu)成線性方程的方法,相應(yīng)地,便確定了每次選擇主元的原則。

然而,在XTX接近奇異時(shí),即det(XTX)近似于0,此時(shí)不管如何確定選取主元的原則,由Tkk變換的性質(zhì)5可知,必然存在某個(gè)主元接近于0,這時(shí)將產(chǎn)生較大誤差。針對(duì)此情況,嶺估計(jì)適當(dāng)增大設(shè)計(jì)陣所有對(duì)角元素的值,使得原本近似于0的特征根變大;并且由于增加量并不大,矩陣本身未作太大變化,但復(fù)共線性降低,此時(shí)得到的參數(shù)估計(jì)在均方誤差意義下優(yōu)于LS估計(jì)?；诖怂枷?本文利用交叉檢驗(yàn)思想,在每次消去變換前適當(dāng)增大主元的大小使得預(yù)測(cè)效果最優(yōu)。

2新算法

在線性模型(1)式中,設(shè)X、Y已作中心化,假設(shè)已將樣本C=(XY)分成3份,即C1=(X1Y1),C2=(X2Y2),C3=(X3Y3)。其中,C1為n1×(p+1)矩陣;C2為n2×(p+1)矩陣;C3為n3×(p+1)矩陣。

2.1數(shù)據(jù)處理

(5)

其中,R(0)表示未作高斯消去變換時(shí)相關(guān)陣的增廣矩陣,而下述出現(xiàn)的R(l)則表示已作l次高斯消去變換后的增廣矩陣。

2.2循環(huán)過(guò)程

(1)主元選擇。計(jì)算所有未選的變量因子的方差貢獻(xiàn),即

(6)

其中,Q={1,2,…,p};L為已選擇添加的變量因子的下標(biāo)集;l為已作高斯消去變換的次數(shù)。顯然,第1次選擇添加變量因子時(shí),L=?,l=0。

(7)

(8)

2.3循環(huán)結(jié)束

(9)

(10)

3算例

現(xiàn)有一樣本數(shù)據(jù)矩陣,將其分成訓(xùn)練樣本和檢驗(yàn)樣本2份。假設(shè)已中心化樣本矩陣,訓(xùn)練樣本C1=(X1Y1),檢驗(yàn)樣本C2=(X2Y2),其中,X1、Y1、X2、Y2分別為:

將訓(xùn)練集樣本轉(zhuǎn)化為相關(guān)矩陣:

此時(shí)得到標(biāo)準(zhǔn)模型回歸系數(shù)為：

4實(shí)驗(yàn)?zāi)M

通過(guò)蒙特卡羅模擬法來(lái)驗(yàn)證新算法的優(yōu)良性,以維數(shù)p=5為例,在總體G:X～Np(0,Σ)中隨機(jī)抽取樣本容量為n的自變量樣本矩陣X,其中協(xié)方差矩陣Σ為:

其中,ρ為自變量樣本的自相關(guān)系數(shù),介于0和1之間,它控制著任意2個(gè)自變量因子的相關(guān)性。

因變量的觀測(cè)值通過(guò)(11)式獲得:

(11)

其中,εi為獨(dú)立正態(tài)分布(0,σ2)偽隨機(jī)數(shù);σ為隨機(jī)誤差的標(biāo)準(zhǔn)差;β=(3.01.50.80.1-2.0)T為固定向量。

考察在σ和ρ的6種不同取值情況下,新算法、LS估計(jì)以及嶺估計(jì)在σ=0.8,1.2,1.5,ρ=0.900,0.999時(shí)的各自估計(jì)效果。

在總體 G中隨機(jī)生成3份樣本容量分別為n1=20,n2=100,n3=200的自變量樣本矩陣X1、X2和X3,按(11)式生成因變量觀測(cè)矩陣Y1、Y2和Y3。

對(duì)得到的3份樣本按上述算法設(shè)計(jì)進(jìn)行計(jì)算。在每個(gè)給定的參數(shù)σ和ρ中,這樣的模擬計(jì)算都重復(fù)進(jìn)行100次。

采用重抽樣法,在3種方法各自得到的100份平均殘差平方和(MRSS)與100份平均預(yù)測(cè)殘差平方和(MPRESS)中,各隨機(jī)抽取500份數(shù)據(jù),并求出中位數(shù)和標(biāo)準(zhǔn)差,見(jiàn)表1和表2所列。

接著生成6種σ和ρ不同取值時(shí)3種估計(jì)法的平均殘差平方和以及平均預(yù)測(cè)殘差平方和的箱線圖,限于篇幅,本文只選取2種有代表性的箱線圖,如圖1、圖2所示。

由表1中可以看到,在重抽樣下,自變量樣本矩陣自相關(guān)性ρ很大(復(fù)共線性很強(qiáng)),新算法的參數(shù)估計(jì)平均殘差平方和的中位數(shù)和標(biāo)準(zhǔn)差都小于嶺估計(jì)。

由表2中可以看到,在重抽樣下,當(dāng)自變量樣本矩陣的自相關(guān)性越強(qiáng),隨機(jī)誤差的標(biāo)準(zhǔn)差σ2越大時(shí),新算法的平均預(yù)測(cè)殘差平方和的中位數(shù)和標(biāo)準(zhǔn)差都要小于LS估計(jì)和嶺估計(jì)。

表1　3種算法模擬數(shù)據(jù)平均殘差平方和(MRSS)的中位數(shù)和標(biāo)準(zhǔn)差比較

注:括號(hào)內(nèi)的數(shù)值為平均殘差平方和的標(biāo)準(zhǔn)差。

表2　3種算法模擬數(shù)據(jù)平均預(yù)測(cè)殘差平方和(MPRESS)的中位數(shù)和標(biāo)準(zhǔn)差比較

注:括號(hào)內(nèi)的數(shù)值為平均預(yù)測(cè)殘差平方和的標(biāo)準(zhǔn)差。

另外從箱線圖中可以看到,在自變量樣本矩陣自相關(guān)系數(shù)足夠大時(shí),新算法估計(jì)的平均殘差平方和以及平均預(yù)測(cè)殘差平方和都小于嶺估計(jì),并且四分位距也小,誤差顯得更為穩(wěn)定。

因此,當(dāng)線性模型的隨機(jī)誤差的方差和自變量樣本自相關(guān)性較大時(shí),新算法下的線性模型參數(shù)估計(jì)的預(yù)測(cè)結(jié)果無(wú)論是精度還是穩(wěn)定性都好于LS估計(jì)和嶺估計(jì),是一個(gè)更為優(yōu)良的估計(jì)。

圖1　σ=1.2,ρ=0.999時(shí)的箱線圖

圖2　σ=1.5,ρ=0.999時(shí)的箱線圖

5結(jié)束語(yǔ)

對(duì)于線性方程組,為了減小最小二乘估計(jì)所得回歸系數(shù)的誤差,需要先做降低病態(tài)系數(shù)矩陣條件數(shù)的預(yù)處理,又由于變量x1,x2,…,xp的量綱可能差異很大,使得XTX不同對(duì)角元素的所需增量差異巨大。本文先將XTX轉(zhuǎn)換成相關(guān)矩陣的形式,利用逐漸添加最關(guān)鍵變量因子的原則選取相應(yīng)的主元,再利用交叉檢驗(yàn)思想,作增大后高斯消去變換；在Matlab中采用蒙特卡洛模擬法進(jìn)行數(shù)值實(shí)驗(yàn)。通過(guò)與LS估計(jì)以及一般嶺估計(jì)比較,新算法在復(fù)共線性較強(qiáng)時(shí)的預(yù)測(cè)效果和穩(wěn)定性都好于另外2種;但交叉檢驗(yàn)的過(guò)程對(duì)樣本數(shù)據(jù)的切分方法存在著人為因素,并且由于過(guò)程相當(dāng)于逐步增加變量建立方程,在運(yùn)算速度上慢于一般嶺估計(jì)。本文只給出一種新算法,嚴(yán)格的理論性質(zhì)將是今后的研究工作方向之一。

[參考文獻(xiàn)]

[1]陳希孺,王松桂.線性模型中的最小二乘法[M].上海:上?？茖W(xué)技術(shù)出版社,2003:1-12.

[2]陳希孺,王松桂.近代回歸分析:原理方法及應(yīng)用[M].合肥:安徽教育出版社,1987:218-226.

[3]張金槐.線性模型參數(shù)估計(jì)及其改進(jìn)[M].第2版.長(zhǎng)沙:國(guó)防科技大學(xué)出版社,1999:68-72.

[4]Liu K.A new class of biased estimate in linear regression[J]. Commun Stat Theory Methods 1993,22:393-402.

[5]Massy W F. Principal components regression in exploratory statistical research[J].J Amer Statist Associ,1965,60: 234-266.

[6]Hoerl A E,Kennard R W.Ridge regression: biased estimation for non-orthogonal problems[J]. Techometrics,1970,12: 55-68.

[7]Al-Hassan Y M.Performance of a new ridge regression estimator[J].Journal of the Association of Arab Universities for Basic and Applied Sciences,2010,9:23-26.

[8]Khalaf G.A proposed ridge parameter to improve the least squares estimator[J].Journal of Modern Applied Statistical Methods,2012,11:443-449.

[9] 殷藝蕓.線性模型中Liu估計(jì)及兩參數(shù)估計(jì)的進(jìn)一步研究[D].重慶:重慶大學(xué),2012.

[10]常新鋒.線性模型參數(shù)有偏估計(jì)的若干研究[D].重慶:重慶大學(xué),2011.

[11]方開(kāi)泰,全輝,陳慶云.實(shí)用回歸分析[M].北京:科學(xué)出版社,1988:70-78.

(責(zé)任編輯張淑艷)

A new algorithm of biased estimation in linear model

LIANG Fei-bao,LIN Tong-hua

(College of Mathematics and Computer Science,Fuzhou University,Fuzhou 350108,China)

Abstract:In order to reduce the prediction residual sum of squares,and based on the method of ridge estimation that increases the diagonal elements of the regression coefficients matrix,a linear model parameter estimation method using the tool of Gaussian elimination transform is proposed. And the data experiment is carried out. Finally,its boxplots about the mean prediction residual sum of squares and the mean residual sum of squares are compared with those of the least squares estimation and the ridge estimation,and the results show that the new algorithm performs better than these two estimations in the estimation precision and stability under certain conditions.

Key words:prediction residual sum of square;parameter estimation;ridge estimation;Gaussian elimination;boxplot

收稿日期：2015-07-22；修回日期：2015-10-13

基金項(xiàng)目：國(guó)家自然科學(xué)基金資助項(xiàng)目(11301084);福建省自然科學(xué)基金資助項(xiàng)目(2014J01010)

作者簡(jiǎn)介：梁飛豹(1963-),男,福建莆田人,福州大學(xué)副教授,碩士生導(dǎo)師.

doi:10.3969/j.issn.1003-5060.2016.06.027

中圖分類號(hào):O212.1

文獻(xiàn)標(biāo)識(shí)碼：A

文章編號(hào)：1003-5060(2016)06-0854-05