徐志昆
(貴州商學院經濟與金融學院,貴陽 550014)
數據缺失現象普遍存在于各應用領域的統計抽樣中。機械原因、人的主觀失誤、歷史局限、有意隱瞞是數據缺失的主導原因。樣本數據的缺失在很大程度上影響分析建模和研究結果的準確性。
長壽命、高可靠性是科技進步帶來的產品發展趨勢,而傳統壽命試驗周期過長,時間成本巨大,加速試驗便成為解決這個問題的極佳選擇。但過去的大部分加速試驗研究都是在沒有數據缺失的前提下進行,若僅因為數據存在缺失就重作試驗是極不現實的,代價非常高昂。
目前針對數據缺失主要有加權法、插補法和構造特殊統計模型三類處理方法。對于單元無回答的情況多采用加權法對缺失值進行補救,而對于項目無回答的處理多采用插補的方法,包括單一插補和多重插補,及根據工程背景構造特殊統計模型。Politz 等(1949)提出了經典的Politz-Simmons調整法[1]。Dempster等[2]首次提出一種使得不完全數據得到有效估計的EM 算法。Rubin 等[3]基于EM 算法,率先提出了多重插補方法。Brick 等[4]提出了最近鄰插補方法,也即樹枝分類的距離函數匹配方法。Liu 等[5]在1994年進一步提出了ECME 等。1998 年金勇進[6]探討了處理缺失數據中對輔助信息的利用問題。2002 年Schafer 等[7]提出的極大似然估計和多重填補法具有較好的處理精度和較廣的應用范圍。王乃生等[8]給出恒定應力加速壽命試驗中數據缺失時的統計方法。2009 年劉寶慧[9]利用回歸插補給出方差分析。楊貴軍等[10]對高相關輔助變量擇優回歸插補法進行了研究。但從20 世紀90 年代初至今,缺乏處理缺失數據的全新思想[11-12]。
針對加速試驗數據一般隨時間呈現單調變化趨勢且精度高的特點,立足于單一插補的角度,提出插值填補法,分別運用Lagrange 插值、三次樣條插值和B 樣條插值,給出缺失數據合理的替補值,達到減小估計量偏差和數據集完整的目的。
三次樣條插值曲線具有良好的性質,在實用中最為普遍。設Δ是[a,b]的一個劃分,則:
若函數S(x)滿足:
(1)S(x) ∈C2[a,b](注:Ck[a,b]表示區間[a,b]上具有k階連續導數的函數集);
(2)S(xi)=f(xi),i= 0,1,…,n;
(3)S(x)在每個子區間[xi,xi+1](i=0,1,…,n- 1)上都是次數不超過三次的多項式,且至少在一個子區間上為三次多項式。則稱S(x)為關于劃分Δ的三次樣條函數。
提出使用三轉角法、三彎矩法和B樣條基函數法完成石英擺片加速退化試驗的缺失值處理。
設S(x)在節點xi(i= 0,1,…,n)處的一階導數值為S'(xi)=mi其中mi是待定參數。記
則有方程組
第二邊界條件S"(a) =f"(a),S"(b) =f"(b),有方程組:
第三邊界條件m0=mn,m1=mn+1,得方程組:
其中:
選擇二階導數作為待定參數:
三彎矩法基本方程[13]:
其中:
在實際應用中,若三次樣條插值沒有邊界條件,最常用的方法就是采用非扭結條件,即:
再由三彎矩基本方程,可得:
利用上述兩種方法均可解出mi(i=0,1,…,n)后,分別代入
即得插值函數S(x),用于插值計算。
B 樣條曲線具有局部性,控制頂點只影響部分曲線的形狀,對其余部分不產生影響,比較具有穩健性,且其造型靈活,還可進行統計數據的光滑化處理。
設有控制頂點P0,P1,…,Pn,則p階(p-1次)B樣條曲線的數學表達式為
其中:Ni,p(μ)是p-1 次B 樣條曲線的基函數。B樣條基函數是一個稱為節點矢量的非遞減的參數μ的序列所決定的p階分段多項式,也即為p階(p-1次)多項式樣條。
B樣條de Boor-Cox遞推定義:
給定一組數據{Qk}(k= 0,1,…,n),找一條p次B 樣條曲線順序通過這組數值點[14],即是根據數據點分布情況選定一組合適的節點矢量U=和控制頂點確定p次B樣條曲線
石英擺片(2010-11-15)在加速應力85℃下試驗的等效撓度數據(單位:10-1mm)如表1 所示,為驗證方法可行性,設定空白處為缺失數據。實際觀測到T1、T2、T6 時刻正面等效撓度分別為7.3856、7.3949、7.3670,T3、T8、T9 時刻反面等效撓度分別為7.3960、7.3900、7.3930。
在實驗中,分別用Lagrange 線性和三次插值、三轉角和三彎矩插值法、均勻和非均勻B樣條基函數法插值進行內推,得到插補結果如表2 所示,并計算誤差平方和(SSE)如表3 和4所示。
插補效果如圖1、圖2所示。
從石英擺片加速退化試驗的兩組數據來看,通過比較圖像和誤差平方和,Lagrange 插值、三轉角法和非均勻B 樣條在兩組缺失數據的插補中均取得了較高的精度,效果理想。

表2 插補結果

圖1 正面等效撓度插補效果

表4 反面等效撓度誤差平方和
石英擺片(2010-01-09)在加速應力85℃下試驗的等效撓度數據(單位:10-1mm)如表5 所示,仍設定空白處為缺失數據。實際觀測到T9、T10 時刻正面等效撓度分別為7.6390、7.6479,反面等效撓度分別為7.6291、7.6359。

表5 石英擺片(2010-01-09)等效撓度
由于插值法一般在外推時精度不高,甚至可能會發生龍格現象,造成巨大偏差。所以在外推缺失數據時,借鑒均值插補方法,將外推轉化為內推處理。步驟如下:
(1)利用已知數據均值來代缺失值相鄰的下一時刻數據;
(2)利用插值法計算缺失數據;
(3)重復步驟(1)、(2)直到所有缺失值計算完成。
得到插補結果如表6所示,并計算誤差平方和(SSE)如表7和表8所示。

表6 插補結果
外推插補效果如圖3、圖4所示。

圖3 正面等效撓度插補效果

圖4 反面等效撓度插補效果

表7 正面等效撓度誤差平方和

表8 反面等效撓度誤差平方和
從加速試驗數據端點缺失的插補結果來看,借助均值插補方法把外推轉化為內推可以避免端點處的巨大波動,降低外推風險,且又在一定程度上反映了數據自身的變化趨勢。通過比較圖像和誤差平方和,Lagrange 插值、三轉角法和非均勻B 樣條在兩組缺失數據的外推插補中得到了較高的精度,效果較好。
數據缺失是統計工作中普遍存在的現象。掌握數據缺失的處理方法,有助于在進行數據采樣、統計分析等環節減少、規避重要信息的丟失,達到提高分析精度的目的。
通過采用Lagrange 插值、樣條插值的單一插補方法來研究加速試驗缺失數據,發現可用于填補缺失數據的中間點,且相比之下算法簡單、容易實現,特別當數據點呈較強規律變化時效果更好。在石英擺片加速退化試驗的缺失數據插補中,插值法從數據點自身變化趨勢出發,并在外推過程中,借鑒均值插值思想把外推轉化為內推,得到了較高精度的缺失數據。
但每一種插值填補方法都不是普遍適用的,都只是對缺失數據分析的一種嘗試。尤其對端點的缺失數據填補應進一步研究。在分析具體問題時,應該綜合權衡考慮使用一種或者幾種方法的綜合結果。