直線擬合的技術(shù)綜合法及穩(wěn)健性判據(jù)

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

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

誤差分析和實(shí)驗(yàn)數(shù)據(jù)處理是物理實(shí)驗(yàn)的重要環(huán)節(jié)。在物理實(shí)驗(yàn)教學(xué)中,當(dāng)已知兩個(gè)物理量成較嚴(yán)格的直線關(guān)系時(shí),常常先測(cè)量多組散布的實(shí)驗(yàn)數(shù)據(jù),再通過(guò)擬合的方法求出直線斜率和截距的最佳估值。直線擬合有多種方法。常用的最小二乘法使殘差平方和極小,計(jì)算簡(jiǎn)便,然而穩(wěn)健性(抗粗差性)較差。最小一乘法使殘差絕對(duì)值之和極小,但算法較復(fù)雜。本文在分析傳統(tǒng)直線擬合方法的特點(diǎn)的基礎(chǔ)上,提出了技術(shù)綜合方法,并通過(guò)蒙特卡羅方法定量地比較了不同直線擬合方法的穩(wěn)健性。

1 傳統(tǒng)直線擬合法的特點(diǎn)分析

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

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

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

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

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

當(dāng)n′為偶數(shù)且n′≤n時(shí),用式由n個(gè)數(shù)據(jù)點(diǎn)(xi,yi)求斜率,是前計(jì)算機(jī)時(shí)代國(guó)外早已使用的簡(jiǎn)化擬合法[5]。當(dāng)高度異常值點(diǎn)出現(xiàn)在直線線段端點(diǎn)附近時(shí),由|?b1/?yi|≈const∝|xi-|0,可看出當(dāng)異常點(diǎn)在直線兩端時(shí)該法穩(wěn)健性可能略優(yōu)于LSM。

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

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

1.5 經(jīng)驗(yàn)調(diào)和法(empirical reconciliation method,EMM)

LSM 求斜率時(shí)有|?b1/?yi|∝|xi-|,SFM求斜率時(shí)有|?b1/?yi|=const,LARM 求斜率時(shí)可粗略地看作各個(gè)yi的貢獻(xiàn)相近。我們?cè)岢鼋?jīng)驗(yàn)調(diào)和法,力圖調(diào)和上述方法的特點(diǎn),對(duì)各yi乘以權(quán)重因子wi=|xi-+ε|-0.5,式中ε是小量。先以為新因變量,以和xi為自變量作“截距”為零的LSM 二元擬合求出b0和b1。再對(duì)擬合結(jié)果中大于的可能異常的值,降低一些權(quán)重，作第2輪加權(quán)擬合。

2 穩(wěn)健擬合的技術(shù)綜合法(technically hybrid method,THM)

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

技術(shù)綜合法的思路與步驟如下:

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

第3步,固定b1,在范圍內(nèi)用最小一乘判據(jù)求出截距b0來(lái)。

3 用蒙特卡羅法判斷穩(wěn)健性

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

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

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

回歸預(yù)報(bào)值的最大誤差的絕對(duì)值|eM|=這是擬合穩(wěn)健性判斷用的參量。

(3) 在l、n與“粗差”點(diǎn)位置相同時(shí),用m組不同的誤差分布數(shù)據(jù)分別計(jì)算標(biāo)準(zhǔn)偏差sy、殘差絕對(duì)值均值與預(yù)報(bào)值最大誤差絕對(duì)值|eM|,再作穩(wěn)健性的統(tǒng)計(jì)分析。

4 用超限率二階矩構(gòu)建擬合穩(wěn)健性的統(tǒng)計(jì)量

4.1 以均勻分布與高斯分布的卷積表征擬合預(yù)報(bào)值的誤差極值分布

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

當(dāng)因變量個(gè)數(shù)n與l=|emax/σ|取某組值時(shí),由蒙特卡羅法計(jì)算出回歸預(yù)報(bào)值的最大誤差絕對(duì)值eM。由于emax可能位于n個(gè)不同點(diǎn)位,一組條件下共有2000n個(gè)eM。

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

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

4.3 對(duì)α=0.01時(shí)emi 分布的誤差限U0.99的討論

對(duì)2000n個(gè)eM遞減排序后,以第20n個(gè)數(shù)作為U0.99,p,以超限率二階矩μ2≈2.17×10-4所對(duì)應(yīng)的eM作為,然后取它們的方均根值分析10×14 組U0.99(n,l)的統(tǒng)計(jì)趨勢(shì),可以發(fā)現(xiàn):

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

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

續(xù)表

4.4 用超限率二階矩構(gòu)建穩(wěn)健性新判據(jù)

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

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

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

4.5 用新判據(jù)比較5種擬合方法的穩(wěn)健性

分別計(jì)算出5種擬合方法、140組n或l不同的pol和μ2之后,統(tǒng)計(jì)計(jì)算結(jié)果的及見(jiàn)表2。

不同方法統(tǒng)計(jì)量的比較以LARM 為對(duì)照,因LARM 是經(jīng)典的穩(wěn)健性方法。140組數(shù)據(jù)匯總結(jié)果表明:THM 的穩(wěn)健性優(yōu)于LARM,它們都顯著優(yōu)于LSM。SFM 的穩(wěn)健性最差。表中EMM 的穩(wěn)健性雖然看起來(lái)略優(yōu)于THM,但是這是以高度統(tǒng)計(jì)離群值的大部分剔除為代價(jià)的。詳細(xì)的計(jì)算表明:當(dāng)離群值位于直線兩端時(shí)THM 的穩(wěn)健型優(yōu)于EMM,THM 保留了LARM 的特征并優(yōu)于后者。

表2 5種擬合法、140組擬合的反預(yù)報(bào)值絕對(duì)值的統(tǒng)計(jì)量與穩(wěn)健性比較

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

EMM 與THM 的穩(wěn)健性優(yōu)于LARM,因?yàn)榍皟煞N方法不同程度地降低了粗差l=|emax/σ|較大的單個(gè)數(shù)據(jù)的權(quán)重,l很大時(shí)使用EMM 擬合相當(dāng)于剔除了高度離群值,THM 在求斜率時(shí)也相當(dāng)于剔除了高度離群值,只在求截距時(shí)用到該值。當(dāng)已知兩變量呈嚴(yán)格直線關(guān)系時(shí),宜用EMM作穩(wěn)健回歸;而當(dāng)變量直線關(guān)系的嚴(yán)格性未知時(shí),宜用THM 穩(wěn)健回歸。

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

圖1 隨l增加而變化

求5種方法的準(zhǔn)偏差sy、平均絕對(duì)殘差分別與LARM 的相應(yīng)量之比,比值的平均值見(jiàn)表4。

表3 固定n 時(shí)的 隨l變化數(shù)表

表3 固定n 時(shí)的 隨l變化數(shù)表

表4 5種方法的準(zhǔn)偏差sy、平均絕對(duì)殘差分別與LARM 的相應(yīng)量之比

表4 5種方法的準(zhǔn)偏差sy、平均絕對(duì)殘差分別與LARM 的相應(yīng)量之比

標(biāo)準(zhǔn)差sy是LSM 極小,THM 與LARM 兩種方法sy的均值之差約0.22%;平均絕對(duì)殘差是LARM 極小。LARM 與THM 兩種方法的均值之差約1.8%?？梢?jiàn)THM 的這兩個(gè)評(píng)價(jià)參量與LARM 相差甚小。THM 保留了LARM 的經(jīng)典統(tǒng)計(jì)的主要特征及優(yōu)點(diǎn)。EMM 的sy及的均值分別位于LSM 與LARM 的相應(yīng)量值之間。

5 技術(shù)綜合法在實(shí)驗(yàn)教學(xué)中的應(yīng)用案例

表5 色散曲線7波長(zhǎng)直線擬合時(shí)不同方法的穩(wěn)健性比較

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

6 結(jié)語(yǔ)

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