穩(wěn)健總體最小二乘法一元線性回歸的相對(duì)有效性探討

劉清，葛永慧

0 引言

線性回歸模型是用于解決數(shù)理統(tǒng)計(jì)問(wèn)題中變量之間關(guān)系的常規(guī)模型，最小二乘法是求解其回歸系數(shù)最經(jīng)典的方法。最小二乘法在求解回歸系數(shù)時(shí)只考慮線性回歸模型中觀測(cè)值含隨機(jī)誤差的情況，總體最小二乘法（TLS）是為了顧及觀測(cè)值和系數(shù)矩陣同時(shí)含有隨機(jī)誤差的狀況。然而無(wú)論是最小二乘法還是總體最小二乘法都不具有抵抗粗差的能力，當(dāng)觀測(cè)值或系數(shù)矩陣中包含粗差時(shí)，參數(shù)估計(jì)的結(jié)果將會(huì)被歪曲[1]。

由于測(cè)量人員的工作疏漏以及儀器精密度不同等各種主觀因素，數(shù)據(jù)采集過(guò)程中粗差會(huì)不可避免地出現(xiàn)。為了剔除或減弱粗差對(duì)參數(shù)估計(jì)的影響，在回歸分析中應(yīng)用穩(wěn)健估計(jì)原理，定義穩(wěn)健回歸的理念。RLS法僅能顧及誤差方程中觀測(cè)向量含粗差的情況，因此，在RLS法的基礎(chǔ)上提出了RTLS法——作為一種可顧及系數(shù)矩陣和觀測(cè)向量中涉及粗差的參數(shù)估計(jì)方法被提出。一些學(xué)者將RTLS法引入到線性回歸模型中，并通過(guò)個(gè)別算例中RTLS法得到比RLS法更小的單位權(quán)中誤差和精度較高的參數(shù)解，從而得出在線性回歸中RTLS法比RLS法更為有效的結(jié)論。

穩(wěn)健最小二乘法（RLS）的研究已取得了大量的研究成果[2]，穩(wěn)健總體最小二乘法（RTLS）的研究也取得了一些成果。楊娟等[3]同樣把基于穩(wěn)健估計(jì)的TLS法應(yīng)用于GPS高程異常擬合問(wèn)題，指出基于穩(wěn)健估計(jì)的TLS法能更好地解決GPS高程擬合模型中控制點(diǎn)已知坐標(biāo)含誤差的問(wèn)題。陳瑋嫻等[4]提出了穩(wěn)健總體最小二乘方法為了解決誤差模型中觀測(cè)值含有粗差的狀況，并以小角度坐標(biāo)轉(zhuǎn)換得出當(dāng)兩套坐標(biāo)均含隨機(jī)誤差且未含粗差時(shí)TLS法優(yōu)于LS法，當(dāng)兩套坐標(biāo)均含隨機(jī)誤差且觀測(cè)值同時(shí)受粗差影響時(shí)RTLS法優(yōu)于TLS法。汪奇生等[5]針對(duì)線性回歸中自變量和因變量可能含有粗差的情況，提出線性回歸模型的穩(wěn)健總體最小二乘算法，并以一元線性回歸為例，得出RTLS法較LS法、TLS法、RLS法能更有效地剔除粗差。

迄今為止，并沒有明確的理論研究說(shuō)明在一元線性回歸中RTLS法的相對(duì)有效性，如果僅僅通過(guò)特殊或極個(gè)別算例就論證在線性回歸中穩(wěn)健總體最小二乘法一定優(yōu)于穩(wěn)健最小二乘法太過(guò)片面，并不足以讓人信服。其次，判斷兩種參數(shù)估計(jì)方法的相對(duì)有效性的比較指標(biāo)應(yīng)選取均方誤差，并不能簡(jiǎn)單地以單位權(quán)中誤差來(lái)衡量精度?；谏鲜龇治觯瑸榱烁_地進(jìn)行測(cè)量數(shù)據(jù)處理和后續(xù)理論的研究指導(dǎo)，對(duì)穩(wěn)健總體最小二乘法在一元線性回歸中的相對(duì)有效性加以分析研究，仍具有一定的迫切性和必要性。

1 模型分類和方法

1.1 線性回歸誤差影響模型的分類

測(cè)量數(shù)據(jù)不可避免地會(huì)受隨機(jī)誤差和粗差的影響，當(dāng)運(yùn)用不同的參數(shù)估計(jì)方法或建立不同的函數(shù)模型時(shí)，使得不同的測(cè)量數(shù)據(jù)在參數(shù)估計(jì)模型中所起的作用也不相同。因此，對(duì)一元線性回歸中觀測(cè)值和系數(shù)矩陣受粗差和隨機(jī)誤差不同影響的變化情形下，劃分為三種不同的誤差影響模型，如下：

模型1：此誤差影響模型中僅觀測(cè)值含有隨機(jī)誤差和粗差，系數(shù)矩陣不含隨機(jī)誤差和粗差。即線性回歸模型中體現(xiàn)為僅因變量含隨機(jī)誤差和粗差。

模型2：此誤差影響模型中僅觀測(cè)值含隨機(jī)誤差，系數(shù)矩陣含隨機(jī)誤差和粗差。即線性回歸模型中體現(xiàn)為自變量含隨機(jī)誤差和粗差，因變量?jī)H含隨機(jī)誤差。

模型3：此誤差影響模型中觀測(cè)值含隨機(jī)誤差和粗差，系數(shù)矩陣僅含有隨機(jī)誤差。即線性回歸模型中體現(xiàn)為因變量含隨機(jī)誤差和粗差，自變量?jī)H含隨機(jī)誤差。

1.2 穩(wěn)健估計(jì)方法

6種常用的穩(wěn)健估計(jì)方法如下[2]：

（1）Huber法：

（2）L1法（殘差絕對(duì)和最小法）：

（4）German-McClure法：

（5）IGG方案：

（6）IGGⅢ方案：

1.3 比較參數(shù)估計(jì)方法的指標(biāo)[2]

（1）殘余真誤差均方誤差（參數(shù)估計(jì)的絕對(duì)指標(biāo)）

其中，設(shè)觀測(cè)值為L(zhǎng)k；觀測(cè)值真值由L?k表示；由參數(shù)估計(jì)方法獲得的Lk的估值由L?k表示；Lk的真誤差由Δk表示；由參數(shù)估計(jì)方法獲得的Lk的改正數(shù)由Vk表示。其

殘余真誤差均方誤差（Mean Square of Residual True Error），由MSRTE表示：

σ?f為殘余真誤差均方誤差，此指標(biāo)能從根本上將參數(shù)估計(jì)方法的優(yōu)劣性和有效性進(jìn)行實(shí)質(zhì)說(shuō)明。用統(tǒng)計(jì)學(xué)的思想對(duì)比兩種參數(shù)估計(jì)方法的相對(duì)優(yōu)劣性和相對(duì)有效性，針對(duì)同一個(gè)參數(shù)估計(jì)問(wèn)題，相同參數(shù)估計(jì)方法進(jìn)行1000次仿真實(shí)驗(yàn)所獲取的MSRTE的平均值定為殘余真誤差均方誤差[6]，由σ?f來(lái)表示。

（2）參數(shù)估計(jì)的MSRTE之比（相對(duì)指標(biāo)）

在同一參數(shù)估計(jì)問(wèn)題中，σ?fa表示由方法A得到的MSRTE，σ?fb表示由方法B得到的MSRTE，令：

R稱為殘余真誤差均方誤差比。式中，是A方法得到的MSRTE，是B方法得到的MSRTE。當(dāng)R＞1.0時(shí)，A方法優(yōu)于B方法；當(dāng)R＜1.0時(shí)，B方法優(yōu)于A方法;當(dāng)R趨向于1.0時(shí)，A和B兩種方法等價(jià)。R作為能從實(shí)質(zhì)上對(duì)兩種參數(shù)估計(jì)有效性根本說(shuō)明的指標(biāo)存在。繼而，運(yùn)用統(tǒng)計(jì)學(xué)的思想對(duì)兩種參數(shù)估計(jì)方法的相對(duì)有效性進(jìn)行說(shuō)明，和通常是指A、B兩種參數(shù)估計(jì)方法對(duì)于同一個(gè)參數(shù)估計(jì)問(wèn)題仿真實(shí)驗(yàn)多次（如1000次）的平均值。

在本文中，RTLS法得到的觀測(cè)值（回歸系數(shù)）估值的殘余真誤差均方誤差與RLS法得到的觀測(cè)值（回歸系數(shù)）估值的MSRTE之比用指標(biāo)RR特指，從根本上說(shuō)明在一元線性回歸中，穩(wěn)健最小二乘法和穩(wěn)健總體最小二乘法的相對(duì)有效性。

2 算例、實(shí)驗(yàn)與討論

2.1 算例

一元線性回歸模型的理論方程為：

從區(qū)間[10,37]中任意選取10個(gè)均勻分布的數(shù)為x，并通過(guò)上述方程獲取真值Y?i，組成10對(duì)點(diǎn)作為理論模擬值。由模擬觀測(cè)真值Y?i、x，在三種不同誤差影響模型下，加上包含粗差ε=10σ0、其余服從正態(tài)分布N(0，0.12)的隨機(jī)誤差得到相應(yīng)的模擬觀測(cè)值yi、xi。

以穩(wěn)健估計(jì)Huber法計(jì)算結(jié)果為例，運(yùn)用RLS法和RTLS法分別進(jìn)行參數(shù)估計(jì)運(yùn)算，得出觀測(cè)值的改正數(shù)V、觀測(cè)值估值的殘余真誤差以及回歸系數(shù)估值a?0、a?1和回歸方程，并通過(guò)觀測(cè)值估值的殘余真誤差fk計(jì)算觀測(cè)值估值的MSRTE以及通過(guò)回歸系數(shù)估值a?0、a?1計(jì)算回歸系數(shù)估值的殘余真誤差。其中，RTLS具體計(jì)算步驟參照文獻(xiàn)[5]。三種誤差影響模型下RLS法和RTLS法的模擬觀測(cè)數(shù)據(jù)和實(shí)驗(yàn)結(jié)果見表1。觀測(cè)值（因變量）真值由Y?i表示，其中元數(shù)由i表示，i=1，2，3。x表示自變量真值，Δ表示服從正態(tài)分布N(0，0.12)的隨機(jī)誤差，y表示由Y?i和隨機(jī)誤差Δ、粗差ε=10σ0共同生成的模擬觀測(cè)值。RLS法的觀測(cè)值改正數(shù)由VA表示，RLS的殘余真誤差由△A表示。RTLS的觀測(cè)值改正數(shù)由VB表示，RTLS的殘余真誤差由△B表示。

表1 模擬觀測(cè)值和真值與RLS法和RTLS法的實(shí)驗(yàn)結(jié)果

由RLS 法計(jì)算得到的回歸系數(shù)估值a?0=3.0233，a?1=0.2481。回歸方程：

由表1中△A列數(shù)據(jù)可計(jì)算得到RLS法觀測(cè)值估值y?的MSRTE為：σ?f1=0.73。

RTLS 法 得 到 的 回 歸 系 數(shù) 估 值a?0=1.0436，a?1=0.3338?；貧w方程：

由表1中△B數(shù)據(jù)得到的RTLS法的y?的MSRTE為：σ?f2=1.05 。

觀測(cè)值估值y?的RTLS法和RLS法的殘余真誤差均方誤差比：

通過(guò)一次計(jì)算就可以獲得觀測(cè)值估值y?的MSRTE，對(duì)于仿真1000次實(shí)驗(yàn)計(jì)算，取其平均值作為觀測(cè)值的MSRTE。

由RLS法計(jì)算得到回歸系數(shù)估值a?0的殘余真誤差為Δ3=a?0-a?0=0.7733 ，回歸系數(shù)估值a?1的殘余真誤差為Δ4=a?1-a?1=-0.0019 。

回歸系數(shù)估值a?0通過(guò)RTLS法得到的殘余真誤差為Δ3=a?0-a?0=-1.2064 ，回歸系數(shù)估值a?1的殘余真誤差為Δ4=a?1-a?1=0.0838 。

回歸系數(shù)估值的殘余真誤差可以通過(guò)一次仿真實(shí)驗(yàn)計(jì)算出，其殘余真誤差均方誤差可由1000次的仿真實(shí)驗(yàn)計(jì)算出。

2.2 一元線性回歸的仿真實(shí)驗(yàn)

當(dāng)給定a?0、a?1和自變量x的取值時(shí)，就可通過(guò)一元線性回歸的理論模型y?=a?0+a?1x得到理論模擬值，通過(guò)運(yùn)用三種不同誤差影響模型，由理論觀測(cè)值加上隨機(jī)誤差或粗差生成模擬觀測(cè)值。

選取三組不同斜率的一元線性回歸方程：A組、B組、C組，分別為斜率約為tan15°、tan45°、tan75°的一元線性回歸模型，相對(duì)應(yīng)的理論回歸方程為：y?=2.25+0.25x、y?=2.25+1.05x、y?=2.25+3.75x。

在上述A、B、C三組方程中分別根據(jù)三種不同誤差影響模型、不同穩(wěn)健估計(jì)方法在觀測(cè)值n=6(n表示觀測(cè)值個(gè)數(shù))的情形下進(jìn)行仿真實(shí)驗(yàn)。

以觀測(cè)值n=6為例，由上述三種不同斜率（tan15°、tan45°和tan75°）的一元線性回歸模型生成的理論模擬值見表2。

表2 不同斜率的理論模擬值（n=6）

理論模擬值加上隨機(jī)誤差或粗差（ε=10σ0）生成模擬觀測(cè)值。對(duì)于觀測(cè)值n=6、三種不同的斜率、三種不同誤差影響模型的情形分別進(jìn)行1000次仿真實(shí)驗(yàn)（方法同算例2.1）。

n=6，ε=10σ0時(shí)RLS和RTLS得到的觀測(cè)值估值y?的殘余真誤差均方誤差見表3，由表3可計(jì)算得RTLS與RLS觀測(cè)值估值y?的殘余真誤差均方誤差之比見表4，回歸系數(shù)的估值a?0和a?1通過(guò)RTLS和RLS所獲得的MSRTE之比見表5。

表3 RLS和RTLS的觀測(cè)值估值y?的殘余真誤差均方誤差(n=6，ε=10σ0)

表4 RTLS和RLS的觀測(cè)值估值 y?的MSRTE之比(n=6，ε=10σ0)

表5 RTLS和RLS的回歸系數(shù)的估值a?0和 a?1MSRTE之比 (n=6，ε=10σ0)

（1）第一種誤差影響模型（僅觀測(cè)值（因變量）含有粗差和隨機(jī)誤差）

對(duì)于A、B、C三組不同斜率（tan15°、tan45°和tan75°）以及觀測(cè)值（因變量）含有粗差和隨機(jī)誤差的情況，通過(guò)RLS和RTLS發(fā)現(xiàn)所得到的觀測(cè)值估值y?的殘余真誤差均方誤差之比（表4）都有RRi＞1.0；RTLS和RLS得到的回歸系數(shù)的估值 a?0和 a?1的殘余真誤差均方誤差之比（表5）是絕大多數(shù)情況下RRi＞1.0和個(gè)別情況下RRi接近1.0。所以RLS無(wú)論是對(duì)于觀測(cè)值估值y?還是對(duì)于回歸系數(shù)的估值a?0和 a?1相對(duì)于RTLS都更有效。

（2）第二種誤差影響模型（系數(shù)矩陣含有隨機(jī)誤差和粗差，觀測(cè)值僅含有隨機(jī)誤差的）

對(duì)于6種常用的相對(duì)有效穩(wěn)健估計(jì)方法，RTLS和RLS在斜率約為tan15°時(shí)，所獲得的觀測(cè)值估值y?的殘余真誤差均方誤差之比（表4）是絕大多數(shù)情況下RRi＞1.0和個(gè)別情況下RRi接近1.0，RLS相對(duì)于RTLS更有效；RTLS和RLS在斜率約為tan45°時(shí)，所獲得的觀測(cè)值估值 y?的殘余真誤差均方誤差之比（表4）1.50＞=RRi＞=0.52，RLS和RTLS有效性很難被說(shuō)明；RTLS和RLS在斜率約tan75°時(shí)，所獲得的觀測(cè)值估值y?的殘余真誤差均方誤差之比（表4）有 RRi＜=0.27；RTLS相對(duì)于RLS更有效。

對(duì)于（tan15°、tan45°和tan75°）這三種不同的斜率以及6種常用的相對(duì)有效的穩(wěn)健估計(jì)方法而言，通過(guò)RTLS和RLS計(jì)算獲取的回歸系數(shù)估值 a?0和 a?1的MSRTE之比（表5）在絕大多數(shù)情況下RRi＞1.0和個(gè)別情況下RRi接近1.0。RLS 對(duì)于回歸系數(shù)的估值 a?0和 a?1比 RTLS 更為有效。

（3）第三種誤差影響模型（觀測(cè)值含有隨機(jī)誤差和粗差，系數(shù)矩陣僅含有隨機(jī)誤差）

對(duì)于A、B、C三組不同斜率（tan15°、tan45°和tan75°）以及6種常用的相對(duì)有效的穩(wěn)健估計(jì)方法而言，RTLS和RLS得到的觀測(cè)值估值 y?的殘余真誤差均方誤差之比（表4）是絕大多數(shù)情形下RRi＞1.0和個(gè)別情況下RRi接近1.0；RTLS和RLS得到的回歸系數(shù)的估值 a?0和a?1的殘余真誤差均方誤差之比（表5）都有RRi＞1.0。即無(wú)論是觀測(cè)值估值 y?還是回歸系數(shù)的估值 a?0和 a?1，RLS相對(duì)于RTLS更有效。

利用類似的方法，對(duì)和n=10和ε=10σ0的情形進(jìn)行了類似的仿真實(shí)驗(yàn)，所得到的結(jié)果與ε=10σ0和n=6的情形相同。

3 結(jié)論

本文用不同的誤差（粗差）影響模型、不同的穩(wěn)健估計(jì)方法和不同的斜率，運(yùn)用仿真實(shí)驗(yàn)的方法，討論了一元線性回歸中RLS法和RTLS法的相對(duì)有效性。

針對(duì)A、B、C三組（斜率約為tan15°、tan45°和tan75°）的一元線性回歸模型：就觀測(cè)值（因變量）的估值而言，在第一種和第三種誤差影響模型下，RLS優(yōu)于RTLS。在第二種誤差影響模型下，當(dāng)斜率較?。stan15°）時(shí)，RLS比RTLS相對(duì)更為有效；當(dāng)斜率約為tan45°時(shí)，難以說(shuō)明RTLS與RLS哪個(gè)更有效；當(dāng)斜率較大（約tan75°）時(shí)，就RLS而言RTLS相對(duì)更為有效。在三種不同的誤差影響模型下，針對(duì)回歸系數(shù)的估值而言，RLS始終都比RTLS相對(duì)更為有效。因此，綜上所述，對(duì)于一元線性回歸RTLS相對(duì)于RLS并沒有明顯的優(yōu)越性。

[1]孫同賀，閆國(guó)慶，周強(qiáng)波.穩(wěn)健初值的選權(quán)迭代法剔除DEM數(shù)據(jù)粗差[J].測(cè)繪科學(xué)，2011,36(3).

[2]葛永慧.再生權(quán)最小二乘法穩(wěn)健估計(jì)[M].北京：科學(xué)出版社，2015.

[3]楊娟，陶葉青.GPS高程異常擬合的文件總體最小二乘算法[J].大地測(cè)量與地球動(dòng)力學(xué)，2014,34(5).

[4]陳瑋嫻，袁慶.抗差總體最小二乘方法[J].大地測(cè)量與地球動(dòng)力學(xué)，2012,32(6).

[5]汪奇生，楊德宏，楊騰飛.線性回歸模型的穩(wěn)健總體最小二乘解算[J].大地測(cè)量與地球動(dòng)力學(xué)，2015,35(2).

[6]葛永慧.再生權(quán)最小二乘法研究[J].測(cè)繪通報(bào),2014,(8).