引入拉格朗日算子的最佳線性回歸模型選擇

閆廣峰 岑敏儀

1 內(nèi)江師范學(xué)院地理與資源科學(xué)學(xué)院，四川省內(nèi)江市東桐路705號，641100 2 西南交通大學(xué)地球科學(xué)與環(huán)境工程學(xué)院，成都市犀安路999號，611756

對線性回歸分析、坐標(biāo)轉(zhuǎn)換和自回歸分析等問題進行建模時，往往需要根據(jù)觀測數(shù)據(jù)的實際情況對模型參數(shù)進行選擇，選取的模型參數(shù)過少會導(dǎo)致所建模型帶有系統(tǒng)誤差，而參數(shù)過多會在一定程度上降低解算精度[1]。由這幾類問題的數(shù)學(xué)模型特點可以發(fā)現(xiàn)，對于平差模型優(yōu)選，均可歸結(jié)為最佳線性回歸模型的選擇問題。對于經(jīng)典線性回歸模型，線性假設(shè)法[2-3]是目前國內(nèi)外公認的模型優(yōu)選有效方法，如利用其確定多項式回歸模型階數(shù)[4]、判斷坐標(biāo)轉(zhuǎn)換模型的相似變換條件是否成立[5]、檢驗三維坐標(biāo)轉(zhuǎn)換中尺度參數(shù)的顯著性[6]以及確定一維自回歸模型和多維自回歸模型階數(shù)[7-8]等。研究表明，利用線性假設(shè)法需滿足一定假設(shè)條件：每次被檢驗的兩個模型中至少存在一個符合觀測數(shù)據(jù)實際，否則得到的結(jié)果不可靠[9-10]。但在實際建模時，每次檢驗的兩個模型是否符合觀測數(shù)據(jù)實際往往未知，因此，采用該方法解決含有多個待選模型的最佳模型選擇問題時，結(jié)果的可靠性和有效性有時無法保證[11]。

此外，經(jīng)典線性回歸模型認為僅因變量包含測量誤差，而忽略自變量的測量誤差。但在實際應(yīng)用中，回歸模型的自變量和因變量均源自觀測數(shù)據(jù)(或測量平差值)，不可避免地含有誤差。因此，傳統(tǒng)回歸分析基于經(jīng)典線性回歸模型并采用最小二乘法求解回歸系數(shù)，其解具有有偏性[12]。對此，同時顧及自變量和因變量測量誤差的回歸模型得到發(fā)展，采用整體最小二乘法求解回歸系數(shù)[13-15]，其解具有漸進無偏性[16]。這些研究成果均以給定的平差模型為基礎(chǔ)，而在實際應(yīng)用中，利用這些方法往往應(yīng)先考慮模型的適用性問題，只有選用合理的平差模型才能夠準確表達觀測量之間的物理或幾何關(guān)系。

目前，除線性假設(shè)檢驗法外，還有白噪聲檢驗準則[10]、殘差平方和準則[1,17]、赤池弘治信息準則[18]等模型優(yōu)選方法，但這些方法同樣僅適用于經(jīng)典線性回歸模型，無法直接用于自變量和因變量均含誤差的回歸模型。由此可見，研究有效的適用于同時考慮自變量和因變量測量誤差的最佳線性回歸模型的選擇理論和方法，是亟待解決的問題。為此，文本基于同時考慮自變量和因變量測量誤差的線性回歸模型，結(jié)合附有參數(shù)約束的測量平差理論和含有多個備選假設(shè)的假設(shè)檢驗理論，研究線性回歸模型的優(yōu)選方法。

1 附有約束條件的線性回歸模型

1.1 顧及自變量和因變量測量誤差的線性回歸模型

線性回歸分析、坐標(biāo)轉(zhuǎn)換和自回歸分析等問題的平差模型均可看作線性回歸模型。在傳統(tǒng)線性回歸分析模型基礎(chǔ)上，為自變量引入誤差向量，可得：

(1)

式(1)中各項誤差的隨機模型為：

(2)

式中，σ0為單位權(quán)中誤差；Qy為觀測量誤差的協(xié)因數(shù)矩陣；設(shè)A中含有u個不同的隨機量，a為將A中不同的隨機量按列依次提取組成的u×1向量，eAs為a對應(yīng)的誤差向量，QAs為a的協(xié)因數(shù)陣。

式(1)可視為非線性高斯-赫爾默特模型。將EA中隨機量作為待估參數(shù)，式(1)可表示為：

(3)

顯然，式(3)為非線性模型，將其轉(zhuǎn)化為線性模型并采用高斯-牛頓法迭代求解。設(shè)第i次迭代后，參數(shù)ξ的估值為ξ(i)，eAs的估值為eAs(i)，將式(3)右端在(ξ(i)，eAs(i))處采用泰勒級數(shù)展開并取至一階項：

(4)

式中，δξ為ξ(i)的微小改正值；R為與ξ(i)及EA有關(guān)的n×u矩陣，滿足ReAs=EAξ(i)。

構(gòu)造目標(biāo)函數(shù)：

(5)

對各變量求偏導(dǎo)并令導(dǎo)數(shù)為0，可得：

(6)

(7)

以參數(shù)的最小二乘解為初值，根據(jù)式(6)和式(7)進行迭代計算，直至‖δξ(i+1)‖2<ε(ε為很小的常數(shù))時迭代結(jié)束。

1.2 附有參數(shù)約束的線性回歸模型

在線性回歸分析、坐標(biāo)轉(zhuǎn)換和自回歸分析等問題對應(yīng)的回歸模型中，參數(shù)個數(shù)少的均可通過對參數(shù)個數(shù)多的附加一定的參數(shù)約束得到。因此，建模時若給定模型參數(shù)選擇的最大范圍，并稱待選模型中參數(shù)個數(shù)最多的為一般回歸模型(general regression model, GRM)，則其他待選模型均可通過對GRM依次添加相應(yīng)參數(shù)約束得到。

假定對某實際問題，共有(f+1)個待選回歸模型，根據(jù)模型中未知參數(shù)個數(shù)由多至少依次排序，設(shè)其分別為F1=0、F2=0…Ff=0、F(f+1)=0。選擇參數(shù)個數(shù)最多的F1=0為GRM，則第j個待選模型Fj=0可通過在式(3)基礎(chǔ)上增加一個約束方程Gjξ=bj得到，即

(8)

式中，j=2,3,…,(f+1)；Gjξ=bj為第j個待選模型對應(yīng)的參數(shù)約束條件。

由式(4)、式(6)和式(7)迭代求解無約束線性回歸模型的系數(shù)矩陣誤差和未知參數(shù)值，設(shè)在第m次迭代后參數(shù)向量滿足收斂條件。

結(jié)合式(4)，將式(8)在(ξ(m)，eAs(m))處采用泰勒級數(shù)展開并取至一階項：

(9)

按求條件極值的方法構(gòu)造目標(biāo)函數(shù)：

(10)

式中，K為對應(yīng)約束條件的聯(lián)系數(shù)向量(拉格朗日算子)。

對各變量求偏導(dǎo)并令導(dǎo)數(shù)為0，可得：

(11)

(12)

拉格朗日算子及其對應(yīng)的方差協(xié)方差陣為：

(13)

(14)

需要說明的是，在實際坐標(biāo)轉(zhuǎn)換問題中，仿射基準變換模型允許不同坐標(biāo)軸采用不同的尺度和旋轉(zhuǎn)參數(shù)，這樣更具一般性，可將其作為GRM；對于線性回歸分析和自回歸分析問題，可首先選擇N/3～2N/3(N為樣本觀測值個數(shù))之間的整數(shù)作為預(yù)設(shè)擬合最高階數(shù)[7]，若在后續(xù)分析中接近預(yù)設(shè)最高階數(shù)仍未得到最佳模型，可再調(diào)整預(yù)設(shè)擬合最高階數(shù)。

2 最佳線性回歸模型選擇

以眾多待選模型中參數(shù)最多的模型作為GRM，則除GRM外的待選模型可表示為附有參數(shù)約束的線性回歸模型，進而可以基于附有參數(shù)約束的測量平差理論和含有多個備選假設(shè)的假設(shè)檢驗理論，從中選出既符合觀測數(shù)據(jù)實際，又具有最少模型參數(shù)的線性回歸模型，稱為初選回歸模型(primary regression model, PRM)。通過分析不同模型平差后的參數(shù)精度，進一步篩選出最佳回歸模型(optimum regression model, ORM)。

2.1 模型初選

要從眾多待選模型中選出最佳模型，首先要對各模型中相應(yīng)的附加參數(shù)約束進行多個備選假設(shè)檢驗和顯著性分析，以定位出不兼容的參數(shù)約束。考慮到最佳模型應(yīng)是參數(shù)解算精度高且模型參數(shù)盡可能少的平差模型[17]，提出原假設(shè)和備選假設(shè)分別為：

(15)

式中，H0為原假設(shè)；H1，H2，…，Hf為f個備選假設(shè)。當(dāng)F(f+1)=0中參數(shù)約束G(f+1)ξ=b(f+1)對回歸模型影響不顯著時接受H0，則F(f+1)=0為PRM；否則拒絕H0，備選模型中PRM需通過進一步假設(shè)檢驗分析來確定。

式(10)中拉格朗日算子是為求解條件極值而引入的中間未知量，已有研究表明[11,19]，其數(shù)值大小反映的是參數(shù)約束對回歸模型影響強度的大小，當(dāng)約束條件對無約束線性回歸模型影響顯著時，拉格朗日算子往往較大，反之較小。基于此，構(gòu)造檢驗統(tǒng)計量：

(16)

(17)

式中，v(m+1)為無約束線性回歸模型在滿足收斂條件時的觀測值殘差向量，Q為其對應(yīng)的協(xié)因數(shù)陣。

給定顯著水平α，其表示H0為真時拒絕H0的概率，只有c′個約束條件全部成立時才接受H0。設(shè)c′個約束條件之間相互獨立，且每個約束條件成立的概率為1-α1，則c′個約束條件全部成立的概率P1為：

(18)

于是有：

1-α=(1-α1)c′=1+c′(-α1)+

(c′(c′-1)/2)(-α1)2+(c′(c′-1)(c′-2)/

(3×2))(-α1)3+…

(19)

由于α1為遠小于1的常數(shù)，因此式(19)可忽略高次項，從而得到：

(20)

α1≈α/c′

(21)

因此對于統(tǒng)計檢驗量T，當(dāng)T

2)將其他待選模型統(tǒng)一為式(8)形式，并在(ξ(m)，eAs(m))處采用泰勒級數(shù)展開，得到式(9)形式。

3)設(shè)原假設(shè)H0和備選假設(shè)HA。H0：待選模型中參數(shù)個數(shù)最少者為PRM；HA：其他模型為PRM。

4)對于由參數(shù)約束最多的模型，根據(jù)式(13)和式(14)得到拉格朗日算子向量K及其對應(yīng)的方差協(xié)方差陣DKK。

5)根據(jù)式(16)和式(21)構(gòu)造假設(shè)檢驗統(tǒng)計量T并進行t檢驗，若T

6)從待選模型組合中刪除原假設(shè)對應(yīng)的線性回歸模型，得到新的待選模型組合，并重復(fù)步驟3)～5)，直到H0成立，算法結(jié)束。

2.2 最佳模型選擇

要從眾多待選模型中得到最佳模型，需要對假設(shè)檢驗得到的PRM與其他待選模型作進一步分析。為此，對比分析所有待選模型平差后的參數(shù)精度。當(dāng)各模型參數(shù)對應(yīng)的中誤差接近時，取參數(shù)個數(shù)最少的模型為最佳模型；而當(dāng)全部或部分模型參數(shù)對應(yīng)的中誤差相差較大時，取最小中誤差對應(yīng)的模型為最佳回歸模型(ORM)。

當(dāng)無約束線性回歸模型取得最優(yōu)解時，待估參數(shù)的方差協(xié)方差陣為：

(22)

當(dāng)?shù)趈個附有參數(shù)約束的待選模型取得最優(yōu)解時，待估參數(shù)的方差協(xié)方差陣為：

(23)

(24)

式中，vj(m+1)為第j個待選模型取得最優(yōu)解時求得的隨機誤差向量，Q為其對應(yīng)的協(xié)因數(shù)陣，c′j為第j個待選模型中參數(shù)約束條件的個數(shù)。

為敘述方便，稱以上包含模型初選和最佳模型選擇等環(huán)節(jié)的線性回歸模型優(yōu)選方法為引入拉格朗日算子的最佳線性回歸模型選擇算法(optimum linear regression model selection algorithm with Lagrange multipliers)，簡稱OLRS-LM法。OLRS-LM算法步驟為：首先將眾多待選模型統(tǒng)一為附有參數(shù)約束的線性回歸模型，然后進行最佳線性回歸模型的篩選分析。由此可見，當(dāng)不考慮回歸模型中自變量的測量誤差時，該方法同樣可行，只是無需進行模型線性化及迭代求解。

3 算例分析

設(shè)計坐標(biāo)轉(zhuǎn)換模型優(yōu)選實驗對OLRS-LM算法的應(yīng)用效果進行分析。選取某GPS網(wǎng)中19個控制點(點號為1～19)在工程獨立坐標(biāo)系下的坐標(biāo)作為Ⅰ套坐標(biāo)系下的模擬真值，并設(shè)計坐標(biāo)轉(zhuǎn)換參數(shù)：平移參數(shù)X0=10 m、Y0=20 m，尺度參數(shù)κx=1、κy=1，旋轉(zhuǎn)角ωx=5.00″、ωy=5.00″，得到控制點在Ⅱ套坐標(biāo)系下的坐標(biāo)真值。控制點兩套坐標(biāo)真值見表1。

在控制點的2套坐標(biāo)中添加期望u=0、中誤差σ=5.0 mm的正態(tài)分布隨機誤差，分別采用OLRS-LM算法和線性假設(shè)法進行坐標(biāo)轉(zhuǎn)換模型優(yōu)選分析。由于經(jīng)典的線性假設(shè)法只能用于傳統(tǒng)線性回歸模型，而且要求被檢驗的2個模型中至少有一個成立，因此，本文在應(yīng)用時對經(jīng)典的線性假設(shè)法進行適當(dāng)改進：1)逐個選擇各待選模型與參數(shù)最多的模型，進行單個備選假設(shè)檢驗；2)對檢驗通過的待選模型兩兩組合進行檢驗，由此得到最佳的平差模型。該過程可確保每次檢驗的2個模型中至少有一個符合觀測數(shù)據(jù)實際，以保證線性假設(shè)法檢驗結(jié)果的可靠性。

在平面坐標(biāo)轉(zhuǎn)換模型中，盡管兩參數(shù)和三參數(shù)模型實際應(yīng)用較少，卻是兩種重要的模型，如在定量分析不同坐標(biāo)系之間的系統(tǒng)誤差時，應(yīng)盡可能地考慮不同坐標(biāo)系之間的差異可能性。因此，將待選轉(zhuǎn)換模型均設(shè)為仿射變換模型、相似變換模型、三參數(shù)模型和兩參數(shù)模型。其中，后3種轉(zhuǎn)換模型均可由仿射變換模型附加相應(yīng)參數(shù)約束得到，具體如下。

1)仿射變換模型：

(25)

式中，(xⅠi,yⅠi)、(xⅡi,yⅡi)分別為第i個點在Ⅰ、Ⅱ兩套坐標(biāo)系中的二維坐標(biāo)，a1=κxcosωx、a2=κysinωy、b1=κxsinωx、b2=κycosωy，(X0，Y0)、(κx，κy)、(ωx，ωy)分別對應(yīng)兩個坐標(biāo)軸的平移、尺度和旋轉(zhuǎn)參數(shù)。

坐標(biāo)變換過程中圖形保持正形特點需滿足柯西-黎曼微分方程：

(26)

式中，Δxi=xⅡi-xⅠi，Δyi=yⅡi-yⅠi。

由式(26)可得滿足正形變換條件的相似變換模型，附加參數(shù)約束為[1]：

(27)

在相似變換模型基礎(chǔ)上，若要得到尺度參數(shù)為1的三參數(shù)模型，附加參數(shù)約束為：

(28)

若2套坐標(biāo)系之間僅有沿坐標(biāo)軸的平移，則為兩參數(shù)模型，附加參數(shù)約束為：

a2=0

(29)

根據(jù)式(25)～(29)可得顧及自變量和因變量測量誤差的仿射變換模型，將其他待選模型表示為附有參數(shù)約束的線性回歸模型。OLRS-LM算法的模型參數(shù)優(yōu)選分析結(jié)果見表2和表3，線性假設(shè)法結(jié)果見表4。其中，2種方法的顯著水平均分別取α1=0.05和α2=0.01。

表4 線性假設(shè)檢驗法最佳模型選擇結(jié)果

從表2和表3可以看出，三參數(shù)模型為初選平差模型，且為最佳平差模型。具體來看，無論取顯著水平α1=0.05還是α2=0.01，對于兩參數(shù)模型，根據(jù)拉格朗日算子構(gòu)造的假設(shè)檢驗統(tǒng)計量遠大于閾值，表明兩參數(shù)模型中附加的部分(或全部)參數(shù)約束與觀測數(shù)據(jù)之間存在不兼容，同時也說明旋轉(zhuǎn)參數(shù)在坐標(biāo)轉(zhuǎn)換模型中具有顯著作用，必須引入；對于三參數(shù)模型，檢驗統(tǒng)計量小于閾值，說明該模型中附加的全部參數(shù)約束與觀測數(shù)據(jù)之間兼容性較好，旋轉(zhuǎn)角在X、Y兩個方向的差異不顯著。分別采用三參數(shù)模型、相似變換模型和仿射變換模型進行平差，兩個平移參數(shù)的中誤差存在明顯差異，其中，三參數(shù)模型和相似變換模型求得的模型參數(shù)精度明顯高于仿射變換模型，且三參數(shù)模型的解算精度高于相似變換模型。由此可以說明，三參數(shù)模型更符合觀測數(shù)據(jù)的實際情況，為最佳平差模型，這與仿真數(shù)據(jù)的實際情況相符。

表2 OLRS-LM算法坐標(biāo)轉(zhuǎn)換模型初選結(jié)果

表3 三種模型平差后坐標(biāo)轉(zhuǎn)換參數(shù)的中誤差

由表4可知，兩參數(shù)模型、三參數(shù)模型和相似變換模型分別與仿射變換模型組合進行差異性檢驗，無論取顯著水平α1=0.05或α2=0.01，兩參數(shù)模型對應(yīng)的檢驗統(tǒng)計量均遠超過閾值，相似變換模型對應(yīng)的檢驗統(tǒng)計量均小于閾值；三參數(shù)模型對應(yīng)的檢驗統(tǒng)計量大于閾值Fα1，小于閾值Fα2。由此說明，Ⅰ、Ⅱ兩坐標(biāo)系之間無顯著的尺度差異，但可能存在旋轉(zhuǎn)角。需要注意的是，此時并無法確定三參數(shù)模型即為參數(shù)解算精度最高的最佳模型，實際上其僅僅是完成最佳平差模型的初選。

分別采用4種待選坐標(biāo)轉(zhuǎn)換模型計算坐標(biāo)轉(zhuǎn)換參數(shù)，并用模型參數(shù)的均方誤差M[21]評定參數(shù)估計結(jié)果的準確度(參數(shù)均方誤差越小，參數(shù)估值準確度越高)，以進一步驗證以上分析結(jié)果的正確性，結(jié)果見表5。

表5 4種待選模型解算得到坐標(biāo)轉(zhuǎn)換參數(shù)及參數(shù)均方誤差

由表5可知，兩參數(shù)模型求得參數(shù)的均方誤差較大，且解算得到的坐標(biāo)轉(zhuǎn)換參數(shù)估計值嚴重偏離設(shè)計真值；而另外3種模型的參數(shù)估計結(jié)果較好，參數(shù)估值與設(shè)計值均相差較小，所得參數(shù)的均方誤差也都較小。三參數(shù)模型、相似變換模型和仿射變換模型求得的坐標(biāo)轉(zhuǎn)換參數(shù)雖然均與設(shè)計值相近，但三參數(shù)模型的均方誤差最小，而且求得的平移、旋轉(zhuǎn)和尺度參數(shù)均與設(shè)計值最接近。由此可見，三參數(shù)模型較相似變換模型和仿射變換模型準確度更高，結(jié)果更優(yōu)。為說明忽略自變量誤差建立回歸模型的不足，進一步建立不考慮自變量誤差的三參數(shù)坐標(biāo)轉(zhuǎn)換模型，并采用最小二乘法求解坐標(biāo)轉(zhuǎn)換參數(shù)，計算其模型參數(shù)均方誤差M，結(jié)果見表6。

表6 忽略自變量誤差的三參數(shù)坐標(biāo)轉(zhuǎn)換模型平差結(jié)果

結(jié)合表5和表6可知，忽略自變量誤差后，所求參數(shù)的均方誤差變大，解算的坐標(biāo)轉(zhuǎn)換參數(shù)除平移參數(shù)Y0基本未變化外，其他轉(zhuǎn)換參數(shù)均與設(shè)計值存在較大差異。

綜上分析可知，只有采用合理的平差模型，才能準確表達觀測量之間的物理或幾何關(guān)系。在實際應(yīng)用中，應(yīng)首先確定最符合觀測數(shù)據(jù)實際的平差模型，然后再進行求解，否則可能得到錯誤的結(jié)果。同時可以發(fā)現(xiàn)，采用OLRS-LM算法和改進的線性假設(shè)法均能對顧及自變量和因變量誤差的線性回歸模型選擇作出客觀、量化的評判，且OLRS-LM算法更為簡便。此外，在相同顯著水平下，采用拉格朗日算子構(gòu)造假設(shè)檢驗統(tǒng)計量可得到正確的判斷結(jié)果，而線性假設(shè)法會出現(xiàn)“棄真”錯誤，表明以拉格朗日算子構(gòu)造假設(shè)檢驗統(tǒng)計量較線性假設(shè)法的檢驗功效更高，得到的結(jié)果更可靠。

4 結(jié) 語

對線性回歸分析、坐標(biāo)轉(zhuǎn)換和自回歸分析問題進行建模時，應(yīng)首先進行模型優(yōu)選分析，以建立既符合觀測數(shù)據(jù)實際、參數(shù)解算精度又高的最佳模型。本文基于附有參數(shù)約束的測量平差理論和含有多個備選假設(shè)的假設(shè)檢驗理論，以拉格朗日算子構(gòu)造假設(shè)檢驗統(tǒng)計量，提出最佳線性回歸模型選擇算法，推導(dǎo)其求解公式，并設(shè)計具體求解算法，通過實驗驗證算法的有效性，得到以下結(jié)論：

1)最佳模型的形式完全由觀測數(shù)據(jù)的實際情況決定，將眾多待選模型統(tǒng)一為附有參數(shù)約束的線性回歸模型，該觀點是線性回歸分析、坐標(biāo)轉(zhuǎn)換和自回歸分析等問題建模時獲得最佳模型的依據(jù)，由此可將模型的優(yōu)選問題轉(zhuǎn)化為含有多個備選假設(shè)的假設(shè)檢驗問題。

2)對于合理的參數(shù)約束，顧及參數(shù)約束的線性回歸模型，其參數(shù)解算精度較無約束的線性回歸模型可得到一定程度提高，因此在實際建模時，對模型進行優(yōu)選分析具有重要意義。

3)OLRS-LM算法可準確找出既符合觀測數(shù)據(jù)實際、參數(shù)解算精度又高的最佳平差模型，其以拉格朗日算子構(gòu)造假設(shè)檢驗統(tǒng)計量，能夠客觀、量化地診斷參數(shù)約束與觀測數(shù)據(jù)之間的兼容性，較以殘差平方和構(gòu)造假設(shè)檢驗統(tǒng)計量的線性假設(shè)法的檢驗功效更高，結(jié)果更可靠。