線性回歸參數(shù)的總體最小二乘估計(jì)新算法

徐 陶，劉國仕，俞 友，吳澤強(qiáng)

（1．湖南省地質(zhì)礦產(chǎn)勘查開發(fā)局402隊(duì)，湖南 長沙410014；2．湖南省勘測設(shè)計(jì)院，湖南 長沙 410014；3．長安大學(xué) 地質(zhì)工程與測繪學(xué)院，陜西 西安710054）

關(guān)于總體最小二乘原理在測量數(shù)據(jù)處理中的應(yīng)用，近年來測量學(xué)者們已進(jìn)行廣泛的研究，并針對(duì)測量平差模型提出一些解決總體最小二乘問題的算法［1－4］，這些算法對(duì)總體最小二乘原理在測量數(shù)據(jù)處理中的應(yīng)用有著重要的作用。但對(duì)于考慮自變量誤差的線性回歸問題，這些算法有一定的局限性。這是因?yàn)榫€性回歸模型中的系數(shù)矩陣存在常數(shù)列，對(duì)其也進(jìn)行改正是不合理的。一般對(duì)于線性回歸的總體最小二乘求解都是采用混合總體最小二乘法［5－7］，但混合總體最小二乘法采用的矩陣分解不利于測量人員理解且沒有考慮到測量平差的優(yōu)勢。文獻(xiàn)［12］先將數(shù)據(jù)中心化，即分離線性回歸系數(shù)矩陣中的常數(shù)列，再采用奇異值分解法進(jìn)行回歸參數(shù)求解。這樣處理比較合理，但計(jì)算較為復(fù)雜。文獻(xiàn)［13］提出一元線性回歸的總體最小二乘迭代算法，但沒有統(tǒng)一到多元線性回歸模型。鑒于此，本文在一元線性回歸總體最小二乘平差模型的基礎(chǔ)上，推導(dǎo)了解線性回歸參數(shù)的總體最小二乘迭代算法，并給出算例，證明算法在進(jìn)行線性回歸參數(shù)估計(jì)時(shí)的正確性和合理性。

1 算法推導(dǎo)

1．1 一元線性回歸總體最小二乘算法

一元線性回歸總體最小二乘平差模型為

式中：y ＝ ［y1… ym］T，x ＝ ［x1… xm］T，vy＝［vy1… vym］T，vx＝ ［vx1… vxm］T。將式（1平差模型展開后按總體最小二乘原理引入平差約束條件為

由式（2）再結(jié)合式（1），則可構(gòu)造目標(biāo)函數(shù)為

由式（3的目標(biāo)函數(shù)可知其含有m＋2個(gè)未知數(shù)，其中2個(gè)為線性回歸參數(shù)a，b，其余的m個(gè)為改正后的自變量＾xi，要求取滿足式（2）的參數(shù)值，即將F分別對(duì)a，b，＾xi求導(dǎo)，并令其等于0。

其中F對(duì)＾xi求導(dǎo)有m個(gè)式子，將其整理化簡后得

F對(duì)a，b求導(dǎo)共有2 m個(gè)式子，將其整理化簡后得

將式（4）寫成矩陣形式為

由于式（5）是F對(duì)回歸參數(shù)a，b求導(dǎo)，而沒有涉及到式（3）中的第二項(xiàng)，因此可以理解為是按最小二乘法求得，只是自變量是改正后的值，而不是原來的觀測值。這時(shí)式（5）可以表示為

式中I＝ ［1 1 … 1］T，為m×1的向量。式（6）、式（7）即是一元線性回歸總體最小二乘迭代算法的基本格式。求解的具體步驟為：①首先按最小二乘原理得到參數(shù)的初值a0，b0；②按式（6）計(jì)算改正后的自變量＾x；③根據(jù)求得的改正后的自變量值由式（7）計(jì)算參數(shù)值；④重復(fù)步驟②和③直到兩次計(jì)算的回歸參數(shù)值之差小于給定的迭代限差，則停止迭代輸出參數(shù)值。

根據(jù)迭代計(jì)算的回歸參數(shù)和自變量改正值代入到式（3）中即可求得殘差平方和，再根據(jù)σ0＝，即可得到單位權(quán)中誤差。

1．2 線性回歸總體最小二乘算法

基于前文的推導(dǎo)，現(xiàn)將一元線性回歸擴(kuò)展到多元線性回歸，其總體最小二乘的平差模型為

式中：

同理將式（8）的平差模型展開后引入總體最小二乘約束條件

由式（9）再結(jié)合式（8）即可構(gòu)造目標(biāo)函數(shù)

從式（10）目標(biāo)函數(shù)中可以看出，其含有m×n＋n＋1個(gè)未知數(shù)。其中的n＋1個(gè)未知數(shù)為回歸參數(shù)a，bj，其余的m×n個(gè)未知數(shù)為改正后的自變量。要求得式（8）中滿足式（9）的一組回歸參數(shù)值，即將F分別對(duì)a，bj，求導(dǎo)并令其等于0。

其中F對(duì)＾xi求導(dǎo)有m×n個(gè)式子，見式（11），將其整理化簡后得式（12）。

F對(duì)a，bj求導(dǎo)共有n＋1個(gè)式子，如式（13），將其整理化簡后得式（14）。

根據(jù)式（11）、式（13）的特點(diǎn)，結(jié)合前節(jié)所述的一元線性回歸總體最小二乘算法推導(dǎo)過程的規(guī)律，可以將其寫成矩陣形式。式（11可表示為

同理式（13）表示為

式中I＝ ［1 1 … 1］T，為m×1的向量，E為n階的單位矩陣。式（15）、式（16）即是多元線性回歸總體最小二乘迭代算法的基本格式。當(dāng)式（8）模型中的n＝1時(shí)即變?yōu)橐辉€性回歸模型，則式（15）、（16）即可變?yōu)槭剑?）和式（7）。因此式（15）、式（16）便是線性回歸總體最小二乘迭代算法的基本格式。

求解的具體步驟為：①按最小二乘原理得到參數(shù)的初值a0，b0；②按式（15）計(jì)算改正后的自變量＾x；③根據(jù)求得的改正后的自變量值由式（16）計(jì)算參數(shù)值；④重復(fù)步驟②和③直到兩次計(jì)算的回歸參數(shù)值之差小于給定的迭代限差，則停止迭代輸出參數(shù)值。

根據(jù)迭代計(jì)算的回歸參數(shù)和自變量改正值代入到式（10）中即可求得殘差平方和，再根據(jù)σ0＝，即可得到單位權(quán)中誤差。

2 實(shí)例分析

為驗(yàn)證本文算法的正確性和可靠性，運(yùn)用Matlab模擬一個(gè)二元線性回歸。其方程為：z＝1．5＋x＋2y，在x和y沒有誤差時(shí)求得z的值上加上均值為0，方差為0．03的隨機(jī)誤差，組成觀測值。然后分別對(duì)x和y添加均值為0，方差為0．03的隨機(jī)誤差，組成新的觀測值，如表1所示。分別采用最小二乘法（LS）、總最小二乘迭代算法、文獻(xiàn)［13］法、本文算法解算線性回歸方程的參數(shù)值，并計(jì)算其單位權(quán)中誤差，結(jié)果如表2所示。

表1 觀測數(shù)據(jù)

續(xù)表1

表2 不同方法解算結(jié)果比較

從表2可以看出，采用總體最小二乘法比采用最小二乘法求得的回歸參數(shù)值更可靠，與真值更為接近，而且精度較高。這是因?yàn)榭傮w最小二乘法考慮了自變量的誤差，提高平差精度。在對(duì)比總體最小二乘迭代算法和本文的總體最小二乘算法時(shí)可以發(fā)現(xiàn)，采用本文算法得到的回歸參數(shù)值與真值最接近，但單位權(quán)中誤差卻較大。這是由于它們對(duì)單位權(quán)中誤差的評(píng)定公式不同。由于常規(guī)總體最小二乘法對(duì)線性回歸模型中的常數(shù)列也進(jìn)行了改正，其改正后的模型為：y＋vy＝（1＋v1）a＋（x＋vx）b，在單位權(quán)中誤差計(jì)算時(shí)將常數(shù)列的改正值當(dāng)作一個(gè)自變量的改正數(shù)，這是不正確的。應(yīng)該將其乘以回歸系數(shù)后的值v1a移項(xiàng)到左邊，作為因變量的另一部分改正數(shù)。如此，計(jì)算的單位權(quán)中誤差為0．028 0，則與實(shí)踐相吻合。另外，采用本文算法得到的結(jié)果與文獻(xiàn)［8］法完全一致，故本文給出的算法正確合理。

3 結(jié)束語

基于線性回歸的總體最小二乘平差模型并以一元線性回歸為基礎(chǔ)，推導(dǎo)了一種迭代算法。該算法的迭代格式簡單，易于編程。相比常規(guī)的總體最小二乘算法，既能考慮到線性回歸模型中系數(shù)矩陣及自變量的誤差，又能顧及系數(shù)矩陣中的常數(shù)列。通過實(shí)例分析，結(jié)果表明針對(duì)線性回歸模型的總體最小二乘問題，本文算法可靠合理。

［1］ GOLUB G H，VAN L C F．An Analysis of the Total Least Squares Problem［J］．SIA M J Nu mer．Anal，1980，17：883－893．

［2］ 魯鐵定，周世健．總體最小二乘的迭代解法［J］．武漢大學(xué)學(xué)報(bào)：信息科學(xué)版，2010，35（11）：1351－1354．

［3］ 許超鈐，姚宜斌，張豹，等．基于整體最小二乘的參數(shù)估計(jì)新方法及精度評(píng)定［J］．測繪通報(bào)，2011（10）：1－4．

［4］ 孔建，姚宜斌，吳寒．整體最小二乘的迭代解法［J］．武漢大學(xué)學(xué)報(bào)：信息科學(xué)版，2010，35（6）：711－714．

［5］ 邱衛(wèi)寧，齊公玉，田豐瑞．整體最小二乘求解線性模型的改進(jìn)算法［J］．武漢大學(xué)學(xué)報(bào)：信息科學(xué)版，2010，35（6）：708－710．

［6］ 丁克良，沈云中，歐吉坤．整體最小二乘法直線擬合［J］．遼寧工程技術(shù)大學(xué)學(xué)報(bào)：自然科學(xué)版 ，2010，29（1）：44－47．

［7］ 孫同賀，羅志才．?dāng)?shù)字化曲線的整體最小二乘平差［J］．工程勘察，2013（8）：71－73．

［8］ 袁豹，岳東杰，張亮，等．基于總體最小二乘的相似變換模型及其在地圖掃描數(shù)字化中的應(yīng)用［J］．測繪工程，2013，22（4）：45－47．

［9］ 李成仁，岳東杰，袁豹，等．基于最小二乘配置的九參數(shù)模型在三維坐標(biāo)轉(zhuǎn)換中的應(yīng)用［J］．測繪與空間地理信息，2014，37（7）：193－196．

［10］龔循強(qiáng)，劉國祥，李志林，等．總體最小二乘擬合問題求解方法的比較研究［J］．測繪科學(xué)，2014，39（9）：29－33．

［11］馮劍橋，黃張?jiān)＃煨憬埽龋傮w最小二乘法在坐標(biāo)轉(zhuǎn)換中的應(yīng)用［J］．測繪與空間地理信息，2014，37（7）：205－206．

［12］邱衛(wèi)寧，陶本藻，姚宜斌，等．測量數(shù)據(jù)處理理論與方法［M］．武漢：武漢大學(xué)出版社，2008．

［13］汪奇生，楊德宏，楊建文．基于總體最小二乘的線性回歸迭代算法［J］．大地測量與地球動(dòng)力學(xué)，2013，33（6）：112－114．