直線擬合的技術綜合法及穩健性判據

郭旭波 肖志剛 朱美紅 常 纓朱鶴年

(清華大學物理系,北京 100084)

誤差分析和實驗數據處理是物理實驗的重要環節。在物理實驗教學中,當已知兩個物理量成較嚴格的直線關系時,常常先測量多組散布的實驗數據,再通過擬合的方法求出直線斜率和截距的最佳估值。直線擬合有多種方法。常用的最小二乘法使殘差平方和極小,計算簡便,然而穩健性(抗粗差性)較差。最小一乘法使殘差絕對值之和極小,但算法較復雜。本文在分析傳統直線擬合方法的特點的基礎上,提出了技術綜合方法,并通過蒙特卡羅方法定量地比較了不同直線擬合方法的穩健性。

1 傳統直線擬合法的特點分析

1.1 最小二乘法(least squares method,LSM)

勒讓德與高斯創立的最小二乘法,使殘差平方和(residual sum of squares,RSS)極小,LSM是統計學發展中的里程碑[1]。直線擬合常用LSM,對n個數據點(xi,yi)求出回歸直線方程=b0+b1xi。LSM 的前提是高斯-馬爾科夫假定:因變量yi具有獨立、同分布的隨機誤差[2,3]。我們可用反證法粗略說明:測量多點散布數據作直線擬合,散布數據擬合主要目的是為減小因變量yi中隨著自變量xi不同而具隨機性的未定系統誤差的影響[4],LSM 令極小而使穩健性(抗粗差性)變差,其中vi=yi-為殘差。直線兩端的數據點yi的誤差ei對擬合斜率b1影響大,因為|?b1/?yi|∝|xi-，其中為xi的平均值。

1.2 最小一乘法(least absolute residuals method,LARM)

穩健性好,它以殘差絕對值之和(sum of absolute residuals,SAR)極小為依據,SAR=由于LARM 算法較復雜,它的出現雖然比LSM 約早40年,但直到60年前進入計算機時代才開始被逐步推廣[1]。

1.3 簡化擬合法(simplified fitting method,SFM)

當n′為偶數且n′≤n時,用式由n個數據點(xi,yi)求斜率,是前計算機時代國外早已使用的簡化擬合法[5]。當高度異常值點出現在直線線段端點附近時,由|?b1/?yi|≈const∝|xi-|0,可看出當異常點在直線兩端時該法穩健性可能略優于LSM。

1.4 先剔除高度異常值后擬合的方法

先剔除高度異常值再回歸是理想的穩健方法,但回歸時的統計離群值判別是尚未系統解決的問題[6,7]。部分文獻中,先用全部n組數據求出回歸直線方程,用最大殘差與其標準差之比G=作為離群值判斷用的統計量。用G＞Gcr,0.99判據、剔除“粗差”后,常出現被剔除值的殘差處于剩余(n-1)組數回歸直線殘差的“包含總體不少于99.9%、概率為99%的統計允許區間”內(依ISO 及GB 標準[8,9]),即sy,n-1,其中sy,n-1為剩余(n-1)組數回歸的標準差,為包含總體不少于99.9%、概率為99%的統計允許區間因子。這說明此類判據往往有邏輯性瑕疵[4]。

1.5 經驗調和法(empirical reconciliation method,EMM)

LSM 求斜率時有|?b1/?yi|∝|xi-|,SFM求斜率時有|?b1/?yi|=const,LARM 求斜率時可粗略地看作各個yi的貢獻相近。我們曾提出經驗調和法,力圖調和上述方法的特點,對各yi乘以權重因子wi=|xi-+ε|-0.5,式中ε是小量。先以為新因變量,以和xi為自變量作“截距”為零的LSM 二元擬合求出b0和b1。再對擬合結果中大于的可能異常的值,降低一些權重，作第2輪加權擬合。

2 穩健擬合的技術綜合法(technically hybrid method,THM)

對于理想模型為截距非零的、自由度不小于5的直線擬合問題,我們提出技術綜合方法,力圖部分綜合上述幾種方法的特點。THM 不作粗差剔除,僅在求斜率時適度降低ki=|yi-大于臨界系數的個別數據yi的權重,以達如下效果:(1)穩健性顯著優于LSM 及SFM,也優于LARM;(2)殘差絕對值和、標準偏差近似與LARM相同;(3)計算過程比LARM 簡單,適用范圍廣。

技術綜合法的思路與步驟如下:

第1步,綜合LSM 與SFM 的特點,用基于LSM 的加權擬合求斜率,使|?b1/?yi|近似正比于|xi-|0.5。

第3步,固定b1,在范圍內用最小一乘判據求出截距b0來。

3 用蒙特卡羅法判斷穩健性

用蒙特卡羅法定量地比較、判斷不同方法的穩健性的步驟為:exp（1）+xiπ的點（xi,yit）。這里設xi等間隔,

(1)先構造n對誤差為零、位于理想直線yit=xi=(2+n)i,i=1,2,…,n。xi在一定區間內呈均勻分布時也能得出與下文類似的結論。

(2) 蒙特卡羅法的參量選擇。建立m組n-1個均值為0、總體標準差σ=0.01的正態分布(或均勻分布)的隨機數組eij(j=1,2,…,m,m=2000),作為n個因變量中n-1 個yi的誤差,yi=yit+eij。剩余一個因變量設置一個人為的“粗差”emax=±lσ,l=3,4,…,16。針對總體正態分布或均勻分布,n取7～16等不同點數,“粗差”點的位置原則上可分別取在直線上的不同位置。分別用LSM、LARM、THM、EMM 及SFM求出直線方程,計算下列特征參量:

回歸預報值的最大誤差的絕對值|eM|=這是擬合穩健性判斷用的參量。

(3) 在l、n與“粗差”點位置相同時,用m組不同的誤差分布數據分別計算標準偏差sy、殘差絕對值均值與預報值最大誤差絕對值|eM|,再作穩健性的統計分析。

4 用超限率二階矩構建擬合穩健性的統計量

4.1 以均勻分布與高斯分布的卷積表征擬合預報值的誤差極值分布

蒙特卡羅法擬合后的每組反預報值的最大誤差為emi,其分布規律不再是正態分布,近似地假設其分布為相同標準差的均勻分布與高斯分布的卷積。設其概率密度為PR&G(em),且有可得99.0%,因此以±2.443σ作為反預報值emi的允許誤差限±LPE,其絕對值LPE≈2.443σ作為|emi|的允許誤差限。

當因變量個數n與l=|emax/σ|取某組值時,由蒙特卡羅法計算出回歸預報值的最大誤差絕對值eM。由于emax可能位于n個不同點位,一組條件下共有2000n個eM。

4.2 超限概率與超限率二階矩

錢鐘泰研究組提出了超限率二階矩(over limited second moment)μ2的概念,其定義式為[4,9]超限率二階矩μ2比超限概率pol能更全面地反映變量e超過誤差限的程度[9],也具有可加性。隨分布類型不同的超限率二階矩μ2有較普遍的可靠性含義。對于分布函數雙側實際上被截尾的情形,μ2反映較多的可靠性信息。誤差呈中心化正態分布時,文獻[9]中給出顯著性水平α=0.01對應的μ2,0.01=2.80×10-4。當emi的分布為相同標準差的均勻分布與高斯分布的卷積時,μ2,0.01≈2.17×10-4,記作。采用±作判斷界限要比±μ2,0.01更苛刻些。

4.3 對α=0.01時emi 分布的誤差限U0.99的討論

對2000n個eM遞減排序后,以第20n個數作為U0.99,p,以超限率二階矩μ2≈2.17×10-4所對應的eM作為,然后取它們的方均根值分析10×14 組U0.99(n,l)的統計趨勢,可以發現:

(2)U0.99能作為判斷穩健性的參量之一,但不同n、l的U0.99(n,l)不具有可加性。如機械地對140組U0.99作平均,結果見表1。從表中可看出,THM 的最小。這與下文及的變化趨勢稍有不同。

表1 5種直線擬合法的U0.99比較

續表

4.4 用超限率二階矩構建穩健性新判據

對2000n個eM遞減排序后,就能依次計算出與eM所對應的超限率二階矩μ2。定出對應誤差限LPE≈2.443σ的μ2,同時定出對應LPE的超限概率pol,它等于eM＞LPE的個數K與2000n之比。μ2等于超出允許誤差限LPE的所有eM的超限率的平方和除以2000n,μ2的含義可以這樣理解:如超限概率為1%,超限這部分eM的超限率約為

比較判斷穩健性程度的統計量要便于對照比較,宜使量綱為1,并在不同參量下具有可加性。我們不用同一置信概率下以某種方法定義的誤差分布的界限值U0.99,而用同一誤差限時的分布特征量μ2。

n、l一定的條件下,當pol＜1%并且μ2＜2.17×10-4時,顯然可認為eM的分布未見異常。對不同n、l的140 組數據匯總后作比較,計算pol＞1%或μ2＞時的總概率以及總超限率二階矩,可以只用作為統計量,以判斷不同擬合方法的穩健性程度。

4.5 用新判據比較5種擬合方法的穩健性

分別計算出5種擬合方法、140組n或l不同的pol和μ2之后,統計計算結果的及見表2。

不同方法統計量的比較以LARM 為對照,因LARM 是經典的穩健性方法。140組數據匯總結果表明:THM 的穩健性優于LARM,它們都顯著優于LSM。SFM 的穩健性最差。表中EMM 的穩健性雖然看起來略優于THM,但是這是以高度統計離群值的大部分剔除為代價的。詳細的計算表明:當離群值位于直線兩端時THM 的穩健型優于EMM,THM 保留了LARM 的特征并優于后者。

表2 5種擬合法、140組擬合的反預報值絕對值的統計量與穩健性比較

對5種方法、140組數的超限概率pol與超限率二階矩μ2值作分析可知:當μ2＞2.17×10-4且pol＞1%時,EMM 只有n=7、l=3 這一組數的μ2明顯大于LSM 的μ2;THM 也只有當n=7、l=3或4這兩組數的μ2明顯大于LSM 的μ2,說明當l＜5時LSM 也有一定的穩健性。EMM所有組的μ2均小于LARM 的μ2;THM 只有當n=9或10、l=5這兩組數的μ2大于LARM 的μ2的（1+1/6）倍,說明EMM 與THM 的穩健性優于LARM。

EMM 與THM 的穩健性優于LARM,因為前兩種方法不同程度地降低了粗差l=|emax/σ|較大的單個數據的權重,l很大時使用EMM 擬合相當于剔除了高度離群值,THM 在求斜率時也相當于剔除了高度離群值,只在求截距時用到該值。當已知兩變量呈嚴格直線關系時,宜用EMM作穩健回歸;而當變量直線關系的嚴格性未知時,宜用THM 穩健回歸。

固定l時,總超限率二階矩隨n的增加大致呈遞減關系。固定n時隨l變化的趨勢見表3。圖1是隨l的增加而變化的曲線。圖中LSM 與LARM 的曲線單調遞增。EMM、THM 與LARM 的曲線大部分在臨界值水平線lg(2.17E4)上方,自下而上依次排列。

圖1 隨l增加而變化

求5種方法的準偏差sy、平均絕對殘差分別與LARM 的相應量之比,比值的平均值見表4。

表3 固定n 時的 隨l變化數表

表3 固定n 時的 隨l變化數表

表4 5種方法的準偏差sy、平均絕對殘差分別與LARM 的相應量之比

表4 5種方法的準偏差sy、平均絕對殘差分別與LARM 的相應量之比

標準差sy是LSM 極小,THM 與LARM 兩種方法sy的均值之差約0.22%;平均絕對殘差是LARM 極小。LARM 與THM 兩種方法的均值之差約1.8%。可見THM 的這兩個評價參量與LARM 相差甚小。THM 保留了LARM 的經典統計的主要特征及優點。EMM 的sy及的均值分別位于LSM 與LARM 的相應量值之間。

5 技術綜合法在實驗教學中的應用案例

表5 色散曲線7波長直線擬合時不同方法的穩健性比較

從表5的最后一行可以看出,即使某一角度錯讀了30′,用穩健性好的技術綜合法擬合,在測量范圍內仍然可得出高置信概率的合理色散曲線,使預報值ni的誤差不大于U0.99。從表5也可看出:雖然THM 的殘差特征值比LARM 大4%或2.4%,但誤差特征值反映的抗“粗差”性較顯著地優于LARM,更是大大優于LSM。

6 結語

本文對直線擬合的技術綜合法的思路和步驟進行了介紹。技術綜合法不作粗差剔除,僅在求斜率時適度降低大于臨界系數的個別數據的權重,且計算過程比最小一乘法簡單,適用范圍廣。用蒙特卡羅法對不同直線擬合方法的穩健性定量比較的結果表明,技術綜合法的穩健性顯著優于最小二乘法,也優于最小一乘法;殘差絕對值和、標準偏差近似與最小一乘法相同。實驗教學應用案例表明,在測三等角棱鏡材料的色散曲線數據時,即使某一角度有較大粗差,用穩健性好的技術綜合法擬合仍然可得出高置信概率的合理色散曲線。