線性回歸參數的總體最小二乘估計新算法

徐 陶，劉國仕，俞 友，吳澤強

（1．湖南省地質礦產勘查開發局402隊，湖南 長沙410014；2．湖南省勘測設計院，湖南 長沙 410014；3．長安大學 地質工程與測繪學院，陜西 西安710054）

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

1 算法推導

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

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

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

由式（2）再結合式（1），則可構造目標函數為

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

其中F對＾xi求導有m個式子，將其整理化簡后得

F對a，b求導共有2 m個式子，將其整理化簡后得

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

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

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

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

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

基于前文的推導，現將一元線性回歸擴展到多元線性回歸，其總體最小二乘的平差模型為

式中：

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

由式（9）再結合式（8）即可構造目標函數

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

其中F對＾xi求導有m×n個式子，見式（11），將其整理化簡后得式（12）。

F對a，bj求導共有n＋1個式子，如式（13），將其整理化簡后得式（14）。

根據式（11）、式（13）的特點，結合前節所述的一元線性回歸總體最小二乘算法推導過程的規律，可以將其寫成矩陣形式。式（11可表示為

同理式（13）表示為

式中I＝ ［1 1 … 1］T，為m×1的向量，E為n階的單位矩陣。式（15）、式（16）即是多元線性回歸總體最小二乘迭代算法的基本格式。當式（8）模型中的n＝1時即變為一元線性回歸模型，則式（15）、（16）即可變為式（6）和式（7）。因此式（15）、式（16）便是線性回歸總體最小二乘迭代算法的基本格式。

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

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

2 實例分析

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

表1 觀測數據

續表1

表2 不同方法解算結果比較

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

3 結束語

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

［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］．武漢大學學報：信息科學版，2010，35（11）：1351－1354．

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

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

［5］ 邱衛寧，齊公玉，田豐瑞．整體最小二乘求解線性模型的改進算法［J］．武漢大學學報：信息科學版，2010，35（6）：708－710．

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

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

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

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

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

［11］馮劍橋，黃張裕，徐秀杰，等．總體最小二乘法在坐標轉換中的應用［J］．測繪與空間地理信息，2014，37（7）：205－206．

［12］邱衛寧，陶本藻，姚宜斌，等．測量數據處理理論與方法［M］．武漢：武漢大學出版社，2008．

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