席夢瑤，賴俊峰，張改梅

（1.內(nèi)蒙古工業(yè)大學(xué)理學(xué)院，呼和浩特 010051；2.呼和浩特市第一醫(yī)院，呼和浩特 010030）

0 引言

在多數(shù)研究中，數(shù)據(jù)缺失現(xiàn)象普遍存在，缺失原因各種各樣。研究數(shù)據(jù)處理不當(dāng)或是無法收集所需信息都會導(dǎo)致數(shù)據(jù)缺失，會對統(tǒng)計分析結(jié)果產(chǎn)生不利影響[1]。為改善該影響，需對原始數(shù)據(jù)進行處理，有多種方法可選擇，如成列刪除、單一插補、多重插補等。且隨著計算能力顯著進步，利用計算機技術(shù)可以產(chǎn)生插補效果更好的方法，如Bootstrap 插補法、機器學(xué)習(xí)插補法等[2]。Bootstrap 插補法在較多領(lǐng)域中均有涉及，熊巍等（2019）[3]指出，針對缺失數(shù)據(jù)可分別采用Bayes法和Bootstrap法進行插補以完成模型構(gòu)建。通過實證對比發(fā)現(xiàn)，兩種方法都有較好的替補效果，Bayes法相較于Bootstrap法結(jié)果更精確，但Bootstrap法更易于操作。Van（2018）[4]針對錯誤分類產(chǎn)生的偏倚，提出用Bootstrap插補法（BI）解決，采用血清肌酐測定來判斷患者是否處于嚴(yán)重腎功能衰竭狀態(tài)，結(jié)果表明BI 避免了無效結(jié)果，在使用準(zhǔn)確疾病概率進行估計時，顯著減小了誤分類偏差。Long 和Johnson（2015）[5]指出，當(dāng)數(shù)據(jù)存在缺失時，變量選擇方法需根據(jù)缺失數(shù)據(jù)機制進行調(diào)整，并提出重抽樣方法（BI-SS），該方法是結(jié)合Bootstrap 插補的方法，大量仿真研究表明，在變量選擇方面，對比其他幾種低維和高維方法，BI-SS 法性能最好，且對調(diào)整參數(shù)值相對不敏感。

概率Weibull 模型最早由數(shù)學(xué)家Weibull提出，用于表示材料斷裂強度的分布，之后用來描述系統(tǒng)或部件具有一定程度可變性的行為。如今，Weibull 分布是描述實際現(xiàn)象概率行為最流行的統(tǒng)計模型之一，在可靠性分析中發(fā)揮著重要作用[6]。對于兩參數(shù)Weibull分布，參數(shù)估計方法主要有極大似然估計法、最小二乘法、矩估計法、對數(shù)矩估計法和L 矩估計法，不同方法條件不同，估計結(jié)果也各不相同[7]。顧蓓青等（2019）[8]研究兩參數(shù)對數(shù)正態(tài)分布和Weibull 分布參數(shù)的近似極大似然估計，通過大量Monte-Carlo 模擬比較這兩種方法估計的精度，結(jié)果表明Weibull 分布參數(shù)的近似極大似然估計結(jié)果更優(yōu)。Tereza和Zuzana（2019）[9]為了比較矩估計法、極大似然估計法、最小二乘法和概率加權(quán)矩估計法對參數(shù)估計的效果，構(gòu)建了57 個模型，并對這些模型用上述方法求解參數(shù)，結(jié)果表明，最小二乘法在44個模型中結(jié)果最優(yōu)，其原因為最小二乘法不需要通過試錯或迭代過程來估計分布的參數(shù)。王晶晶和孫志絢（2022）[10]提出針對小樣本兩參數(shù)Weibull分布的迭代極大似然估計法。根據(jù)數(shù)據(jù)對其進行可靠性建模，使用該方法估計Weibull分布模型參數(shù)，得到可靠性函數(shù)，仿真試驗結(jié)果表明，使用該方法計算的結(jié)果與實際結(jié)果基本一致，具有較高的準(zhǔn)確性與實用性。Liu等（2022）[11]針對滾動軸承壽命的可靠性評估，基于相關(guān)系數(shù)優(yōu)化法和極大似然估計法提出一種估計參數(shù)的方法——Cor-MLE，得到較為準(zhǔn)確的Weibull分布參數(shù)，并用Monte-Carlo模擬試驗驗證該方法的可行性。與最小二乘法相比，該方法具有更高的精度，能滿足滾動軸承的壽命分析。Feng 等（2022）[12]為了推斷產(chǎn)品在正常應(yīng)力條件下的可靠性，假設(shè)數(shù)據(jù)服從Weibull 分布，且假定該分布具有恒定的形狀參數(shù)，尺度參數(shù)則采用極大似然估計法計算，該過程由Newton-Raphson 算法完成，試驗結(jié)果證實了極大似然估計量的有效性，且得到的產(chǎn)品可靠性符合真實情況。Freels等（2019）[13]指出，Weibull分布一直是模擬各種機械和生物壽命數(shù)據(jù)的常用方法，在某些情況下，初始值非常影響迭代收斂性，需對Weibull分布進行修正，可以使用不同方法修正Weibull 分布，用極大似然估計法求解參數(shù)并對比概率密度函數(shù)，結(jié)果表明，修正后的Weibull 分布比未修正的Weibull分布能更好地擬合數(shù)據(jù)集。

產(chǎn)品質(zhì)量和生產(chǎn)效率是企業(yè)能夠穩(wěn)定生存的重要法寶，在保證質(zhì)量的前提下提高效率是企業(yè)的不斷追求[14]。織造是制造織物通用的方法，隨著科學(xué)技術(shù)的發(fā)展，織造速度大幅提高，織物原材料的質(zhì)量帶給高技術(shù)產(chǎn)品的瓶頸效應(yīng)日益突出[15]。經(jīng)紗作為織物的重要原材料之一，影響其質(zhì)量的因素主要有材料、張力等。評價經(jīng)紗質(zhì)量的一個指標(biāo)是經(jīng)紗上機時的斷裂次數(shù)。斷裂是在織造過程中由薄弱環(huán)節(jié)引起的。若要提高生產(chǎn)效率，解決方法之一是減少設(shè)備停機時間，從而減少經(jīng)紗斷裂次數(shù)。如果能挑選出斷裂次數(shù)最少的經(jīng)紗，就可以提高勞動生產(chǎn)率，進而提高織物質(zhì)量[16]。

為厘清采用Bootstrap 插補法處理缺失值是否影響含缺失數(shù)據(jù)Weibull 分布的參數(shù)估計，本文首先對含缺失數(shù)據(jù)Weibull 分布的參數(shù)進行極大似然估計和最小二乘估計，然后采用MCMC 法對Weibull 分布數(shù)據(jù)進行模擬試驗和經(jīng)紗斷裂數(shù)據(jù)驗證，最后根據(jù)試驗結(jié)果進行總結(jié)與展望。

1 含缺失數(shù)據(jù)Weibull分布的參數(shù)估計

1.1 Weibull分布數(shù)據(jù)的Bootstrap處理

假設(shè)一組數(shù)據(jù)Y=(Y1，Y2，…，Yp)，變量Yj表示為Yj=，n=1，2，…，p表示變量個數(shù)，將Y稱為資料陣，用矩陣表示為：

假定a為變量Yj是否缺失的指示變量，當(dāng)a=1 時表示項目有回答，當(dāng)a=0 時表示項目無回答。A為缺失數(shù)據(jù)指示變量數(shù)據(jù)集，用矩陣表示為：

缺失數(shù)據(jù)模式表達了在數(shù)據(jù)集中哪些數(shù)據(jù)被觀測到，哪些數(shù)據(jù)缺失，有助于認(rèn)識數(shù)據(jù)集中不同變量之間的相關(guān)關(guān)系，幫助目標(biāo)變量選擇合適的數(shù)據(jù)處理方法[17]。缺失數(shù)據(jù)模式共分四種，如圖1所示，分別為單變量缺失、多變量缺失、單調(diào)缺失和一般缺失。

圖1 缺失數(shù)據(jù)模式

缺失數(shù)據(jù)機制表達了缺失數(shù)據(jù)與數(shù)據(jù)集變量值之間的關(guān)系。缺失數(shù)據(jù)機制從本質(zhì)上說明了數(shù)據(jù)是如何缺失的，隨后的模型建立都是基于缺失數(shù)據(jù)機制的假定[18]。一組數(shù)據(jù)缺失主要分為三種機制：

（1）完全隨機缺失（MCAR）

其中，L表示分布。

（2）隨機缺失（MAR）

其中，Yobs表示只含觀測數(shù)據(jù)。

（3）非隨機缺失（NMAR）

變量Yj的缺失與自身有關(guān)，不管其是否觀測到，都為不可忽略的缺失。

假設(shè)一組數(shù)據(jù)X=(x1，x2，…，xn) 服從Weibull 分布W(η，β)，其中η和β分別為其尺度參數(shù)和形狀參數(shù)，其分布函數(shù)為：

概率密度函數(shù)為：

可靠性函數(shù)為：

則效率為：

當(dāng)β＞1時，r(x)遞增；當(dāng)0 ＜β＜1時，r(x)遞減[9]。

若數(shù)據(jù)X存在缺失值，按有無缺失值排序，前n1個觀測值未缺失，記為Xobs,后n0個觀測值缺失，記為Xmis，即X=(Xobs，Xmis),且n1+n0=n。Bootstrap 插補法的基本思想是根據(jù)觀測部分Xobs可放回地抽取樣本量為h的數(shù)據(jù)，記為，對其求均值μh，即：

將μh作為后n0個缺失值的插補值，則原始缺失數(shù)據(jù)變?yōu)閄B=(Xobs，μh)，此時均值和方差按正常無缺失值計算，分別為：

由于Bootstrap樣本只是原樣本的一個重復(fù)隨機抽樣，因此沒有損失太多原始信息。Bootstrap 樣本可以無限抽取，且彼此之間存在差異，此差異恰好反映了缺失數(shù)據(jù)的隨機性，彌補了其他方法對總體方差的低估[19]。本文選擇均值μh作為插補值，也可選擇中位數(shù)或眾數(shù)作為插補植，需具體問題具體分析。

1.2 Weibull分布參數(shù)的極大似然估計

極大似然估計是一種用于參數(shù)估計的統(tǒng)計學(xué)方法，由于比其他估計方法簡單，且計算復(fù)雜程度較低，因此本文選擇極大似然估計法來估計Weibull分布的參數(shù)。具體步驟為：

（1）寫出模型的對數(shù)似然函數(shù)；

（2）求對數(shù)似然函數(shù)對各個參數(shù)的偏導(dǎo)，令偏導(dǎo)等于0；

（3）解方程組求參數(shù)的極大似然估計值[20]。

具體公式如下：

Weibull分布似然函數(shù)的表達式為：

整理有：

對數(shù)似然函數(shù)為：

對式（11）分別求關(guān)于η和β的偏導(dǎo)，令其等于0，即：

解方程組（12）可得到η和β的估計值η?和β?，但該方程組為超越方程組，不易求出解析解，故用數(shù)值解替代，選用Newton-Raphson算法求其近似值。

根據(jù)Newton-Raphson 算法的基本步驟，在兩參數(shù)Weibull分布中，先確定非線性方程組：

再確定對應(yīng)的Jacobi矩陣，設(shè)θ=(η，β)，有：

其中：

若矩陣F'(θ)非奇異，則可解出θ的近似值，并把它作為下一次迭代值。于是，Newton-Raphson 迭代公式為：

其中，θ(0)是給定的初始值。在選取一定范圍的精度和初始值后，可完成對方程組（13）的求解[21]，進而得到η?和β?的值。

當(dāng)Jacobi 矩陣非奇異時，可近似得到參數(shù)的估計值，但如果Jacobi 矩陣奇異，或初值不滿足迭代條件，那么結(jié)果將不收斂，無法得到參數(shù)估計值的近似值。此時選擇最小二乘法求解參數(shù)η和β的估計值。

1.3 Weibull分布參數(shù)的最小二乘估計

對可靠性函數(shù)（見式（5））取對數(shù)，有：

對式（15）再取對數(shù)，可得：

令：

則式（16）變?yōu)椋?/p>

中點估計為：

從而得到R(xi)的值，由式（5）可以得到一組數(shù)組{(x1，，經(jīng)過方程組（17）的變形，數(shù)組變?yōu)椤?/p>

定義離差平方和為：

所謂最小二乘估計，就是尋找參數(shù)η和β的估計值，使離差平方和達到最小。對Q分別求關(guān)于參數(shù)η和β的偏導(dǎo)，并令其等于0，從而得到參數(shù)的估計為：

2 模擬驗證

2.1 模擬數(shù)據(jù)

為驗證Bootstrap插補法是否對Weibull參數(shù)估計值有影響，模擬一組尺度參數(shù)為1、形狀參數(shù)為2、樣本量為100的Weibull分布W（1，2），組成完整數(shù)據(jù)集。模擬數(shù)據(jù)隨機缺失10%、25%、40%和60%，隨機缺失不同比例可根據(jù)二項分布實現(xiàn)，分別對缺失數(shù)據(jù)采用Bootstrap插補法進行處理，設(shè)抽樣量分別為10、50、100、500 和1000，計算不同條件下的參數(shù)估計值并進行對比，并且用MCMC法模擬驗證不同缺失比例下Bootstrap插補法的性能。

數(shù)據(jù)隨機缺失不同比例時，缺失情況如圖2所示。其中，淺色代表數(shù)據(jù)觀測值，深色代表數(shù)據(jù)缺失值，對比發(fā)現(xiàn)，隨著缺失比例的增加，缺失數(shù)目逐漸增多，當(dāng)缺失比例為60%時，缺失數(shù)目達到70個。

圖2 不同缺失比例下的數(shù)據(jù)缺失情況

不同情況的缺失數(shù)據(jù)用不同樣本量的Bootstrap插補，采用Newton-Raphson 算法計算參數(shù)估計值的近似數(shù)值解，將采用最小二乘法計算的參數(shù)估計值作為初始值開始迭代，得到對比結(jié)果如表1所示。

表1 不同缺失比例下不同樣本量Bootstrap插補的參數(shù)估計

表1 表明，在不同缺失比例下，參數(shù)估計值有較大差別，但抽樣量對參數(shù)估計值影響較小，主要為缺失比例對參數(shù)估計值的影響。由于當(dāng)缺失比例為60%時，用Newton-Raphson 算法得到的結(jié)果并不收斂，因此選用最小二乘法計算參數(shù)估計值。除缺失比例為60%的情況外，其余缺失比例下的參數(shù)估計值結(jié)果均在允許的誤差范圍內(nèi)，當(dāng)缺失比例為25%時，結(jié)果最接近真實參數(shù)值。

獲得參數(shù)估計的近似值后，需對比不同缺失比例下Bootstrap 插補法的性能，即觀察參數(shù)的偏差和均方誤差。由于抽樣量對參數(shù)估計值的影響較小，因此只對比抽樣量為1000，缺失比例分別為10%、25%、40%和60%的Bootstrap處理后的“完整”數(shù)據(jù)。對每個“完整”數(shù)據(jù)采用MCMC 法計算偏差和均方誤差，且與未缺失數(shù)據(jù)計算結(jié)果對比，如果前者結(jié)果和后者差別不大，那么可以認(rèn)定Bootstrap插補法可作為處理Weibull分布缺失數(shù)據(jù)的方法。

對比參數(shù)β時，控制參數(shù)η不變，變?yōu)閱螀?shù)Weibull分布。偏差bn定義為估計量與真實值之差，即：

一般地，參數(shù)估計量是不唯一的。對同一個參數(shù)的不同估計量，哪一個較為有效，通常情況是利用估計量的均方誤差來判斷。均方誤差是估計量與真實值差值平方的期望，即均方誤差為：

均方誤差越小，有效性越高；均方誤差越大，有效性越低[23]。為計算偏差，先用MCMC 法模擬出未知的。具體操作如下：

（1）將不同缺失比例下Newton-Raphson 算法的參數(shù)估計值作為Weibull 分布的β，模擬樣本量為100 的Weibull分布W；

（3）根據(jù)式（23）和式（24）計算bn和MSE，參數(shù)η的bn和MSE也通過類似方法得到。

根據(jù)上述步驟，依次計算參數(shù)β和η的偏差和均方誤差，計算結(jié)果如表2、圖3和下頁圖4所示。

表2 不同缺失比例下參數(shù)估計的bn 和MSE

圖3 不同情況下用MCMC法計算β 的bn 和MSE

圖4 不同情況下用MCMC法計算η 的bn 和MSE

對于β來說，隨著缺失比例的增加，bn和MSE的值逐漸增大，未缺失數(shù)據(jù)計算的bn為0.0197，MSE為0.0212，差別最大的結(jié)果是缺失比例為60%的“完整”數(shù)據(jù)集，bn和MSE分別為0.0321和0.0566，但總體對比看，和未缺失的計算結(jié)果只相差0.0124和0.0354，屬于合理的誤差范圍。當(dāng)缺失比例為25%時，計算結(jié)果小于缺失比例為10%的結(jié)果，由于數(shù)據(jù)插補具有隨機性，因此以上結(jié)果可以說明本次插補值更接近原始數(shù)據(jù)。對于η來說，整體也是隨著缺失比例的增加，bn和MSE的值逐漸增大，除了缺失比例為60%，在其他缺失比例下，數(shù)值結(jié)果都是等于或小于未缺失數(shù)據(jù)計算的結(jié)果，說明Bootstrap插補法對Weibull分布缺失數(shù)據(jù)的處理是有效的。

2.2 經(jīng)紗斷裂數(shù)據(jù)

當(dāng)經(jīng)紗斷裂數(shù)據(jù)在機器上未能完整測量時，數(shù)據(jù)將缺失一部分，若直接進行統(tǒng)計分析，將導(dǎo)致結(jié)果有誤差。斷裂數(shù)據(jù)通常服從Weibull 分布，引用R 內(nèi)置數(shù)據(jù)集warpbreaks，變量為breaks、tension和wool，表示2種材料的經(jīng)紗在3種不同張力的條件下在機器上運作時的斷裂次數(shù)，一個機器記錄9次，共有54行數(shù)據(jù)。根據(jù)斷裂次數(shù)來擬合曲線，結(jié)果如圖5 所示。由圖5 可知，曲線大致符合Weibull分布曲線，因此假定該組斷裂數(shù)據(jù)服從兩參數(shù)Weibull 分布，記為W(η，β)。為驗證上述試驗的正確性，考慮到斷裂次數(shù)本身并不含缺失值，因此讓數(shù)據(jù)隨機模擬缺失10%和25%，分別對缺失數(shù)據(jù)采用抽樣量為100 和1000 的Bootstrap 插補法處理，由于斷裂數(shù)據(jù)方差過大，為174.204，若對處理后的“完整”數(shù)據(jù)采用極大似然估計法來估計參數(shù)η和β，Newton-Raphson 算法的結(jié)果將不收斂，因此采用最小二乘法求解。參數(shù)η和β的估計結(jié)果如表3所示。

表3 不同缺失條件下斷裂次數(shù)的參數(shù)估計

圖5 斷裂次數(shù)擬合曲線

根據(jù)表3觀察到，抽樣量對插補后的參數(shù)估計值影響較小，影響較大的是缺失比例。當(dāng)缺失比例為10%時，β的估計值與真實值相差0.346，小于0.5，在可接受范圍內(nèi)，但η的差值絕對值達到1.777，是誤差較大的結(jié)果。且缺失比例越高，β的差值絕對值也越大，但Bootstrap 插補具有隨機性，每次插補值不確定，故當(dāng)缺失比例為25%時，η的估計值與真實值相差1.146（n=100）和1.137（n=1000），小于缺失10%的情況?？蛇M一步對比概率密度曲線，由于抽樣量影響不大，因此只對比不同缺失比例下抽樣量為1000 的概率密度曲線，如圖6 所示，可觀察到當(dāng)缺失比例為10%時，其概率密度曲線和真實數(shù)據(jù)相差不大，但當(dāng)缺失比例為25%時，概率密度曲線波動較大，與真實數(shù)據(jù)相差較大，與上述模擬試驗結(jié)果基本吻合。

圖6 當(dāng)數(shù)據(jù)分別缺失10%和25%時,與實際數(shù)據(jù)對比的概率密度曲線

3 總結(jié)與展望

本文探討了服從Weibull分布的數(shù)據(jù)在不同缺失比例下用Bootstrap插補法處理缺失數(shù)據(jù)的效果，用極大似然估計法和最小二乘法求解兩參數(shù)Weibull 分布的估計值，結(jié)果表明兩種方法都可以得到較好的參數(shù)估計值，但用Newton-Raphson 算法計算比最小二乘法誤差小，而Newton-Raphson 算法對數(shù)據(jù)要求較高，故在條件滿足時選擇用Newton-Raphson 算法，條件不滿足時選擇用最小二乘法求解參數(shù)近似解。Bootstrap 插補法可作為處理Weibull分布缺失數(shù)據(jù)的一種方法，由于Bootstrap插補法只是采用對數(shù)據(jù)集的放回抽樣，因此避免了不合理的值。但隨著缺失比例的增加，參數(shù)估計值與真實值的差別逐漸擴大。在大樣本的情況下，且缺失比例不超過60%時，參數(shù)估計值都為合理范圍值，而且真實數(shù)據(jù)缺失一般都不會超過60%，故Bootstrap 插補法可作為處理Weibull 分布缺失數(shù)據(jù)的方法。

調(diào)查中數(shù)據(jù)缺失的變量往往有多個，故插補方法應(yīng)用廣泛。根據(jù)數(shù)據(jù)特點選擇合適的插補方法對整個試驗過程有事半功倍的效果。但本文只進行數(shù)據(jù)模擬研究，且實際經(jīng)紗數(shù)據(jù)距現(xiàn)在時間較長，下一步希望找到合適的調(diào)查數(shù)據(jù)，讓Bootstrap插補法在實際中得到應(yīng)用。