模型估計(jì)對序數(shù)響應(yīng)輪廓控制圖的影響

王志瓊 李金夢 王 穎 馬彥輝

天津理工大學(xué)管理學(xué)院，天津，300384

0 引言

統(tǒng)計(jì)過程控制(statistical process control，SPC)是一種將統(tǒng)計(jì)技術(shù)應(yīng)用于過程監(jiān)控的方法，該方法通過控制圖等工具監(jiān)控過程的穩(wěn)定性，以確保最終產(chǎn)品的質(zhì)量[1]。在實(shí)際中，某些產(chǎn)品或過程的質(zhì)量特性需要用函數(shù)關(guān)系表示，如產(chǎn)品故障率與使用時間的關(guān)系，這種響應(yīng)變量與一個或多個解釋變量之間的函數(shù)關(guān)系稱為輪廓(profile)[2]。對由一系列輪廓組成的時間序列數(shù)據(jù)進(jìn)行監(jiān)控被稱為輪廓監(jiān)控[3]。根據(jù)輪廓響應(yīng)變量的性質(zhì)，輪廓可分為兩種類型：一種是具有數(shù)值響應(yīng)變量的輪廓，另一種是響應(yīng)變量為屬性(計(jì)數(shù)或分類)數(shù)據(jù)的輪廓。對于第二種輪廓，大部分監(jiān)控方法沒有考慮響應(yīng)屬性之間的序數(shù)信息(例如，好、中和差)，而這些信息對構(gòu)建控制方案具有重要的作用。NOOROSSANA等[4]通過使用序數(shù)邏輯回歸擬合輪廓模型，提出了4種方法來監(jiān)控序數(shù)輪廓，這是對序數(shù)輪廓的首次研究。HAKIMI等[5]使用序數(shù)列聯(lián)表擬合序數(shù)變量之間的關(guān)系，并提出了一種基于序數(shù)正態(tài)統(tǒng)計(jì)量的控制圖，用于監(jiān)控序數(shù)對數(shù)線性模型。DING等[6]假設(shè)存在一個潛在的連續(xù)變量決定響應(yīng)變量的水平，提出了一種新的控制方案，用于監(jiān)控參數(shù)回歸模型。本文針對響應(yīng)變量為序數(shù)數(shù)據(jù)的輪廓，采用非參數(shù)模型擬合函數(shù)關(guān)系，并使用經(jīng)典的廣義似然比(generalized likelihood ratio，GLR)控制圖對該函數(shù)關(guān)系進(jìn)行監(jiān)控。

輪廓監(jiān)控一般分兩個階段進(jìn)行，即階段Ⅰ和階段Ⅱ。階段Ⅰ是以采集到的歷史樣本數(shù)據(jù)為基礎(chǔ)，剔除失控(out-of-control，OC)樣本點(diǎn)，以確保剩下樣本均處于受控(in-control，IC)狀態(tài)，從而建立穩(wěn)定的IC模型。階段Ⅱ是在階段Ⅰ建立的IC模型的基礎(chǔ)上對實(shí)時數(shù)據(jù)進(jìn)行監(jiān)控，以及時發(fā)現(xiàn)過程異常[1]。階段Ⅱ控制圖的構(gòu)建往往是假設(shè)IC模型已知或者可以被完全精確估計(jì)[7-8]。然而，在實(shí)際中，模型通常是未知的，需要通過階段Ⅰ的分析，剔除異常點(diǎn)和變點(diǎn)以確保樣本處于IC狀態(tài)，然后進(jìn)行模型估計(jì)。因?yàn)槟Ｐ凸烙?jì)值具有波動性，所以在監(jiān)控過程中使用估計(jì)值代替已知模型可能會影響階段Ⅱ控制圖的性能。在根據(jù)階段Ⅰ的IC輪廓數(shù)據(jù)集估計(jì)模型時，MAHMOUD[9]評估了監(jiān)控簡單線性輪廓的方案的性能。YAZDI等[10]比較參數(shù)估計(jì)對監(jiān)控多元簡單線性輪廓的3種階段Ⅱ方法性能的影響。更多詳細(xì)內(nèi)容見文獻(xiàn)[11-17]。本文的主要目的是探討在IC模型未知且被估計(jì)量取代的情況下，模型估計(jì)對階段Ⅱ序數(shù)輪廓控制圖性能的影響。

首先，樣本量在模型估計(jì)中起較大的作用， JONES等[18]討論了樣本量大小對控制圖性能的影響，認(rèn)為要獲得與模型已知時類似的結(jié)果，需要較大的樣本量。其次，可以選用不同的估計(jì)方法，主要分為兩大類：參數(shù)方法[6]和非參數(shù)方法[19]。本文主要使用兩種非參數(shù)方法即局部線性核估計(jì)(local linear kernel estimation，LLKE)和樣條法(spline)，以及參數(shù)方法Newton Raphson。最后，在不同的估計(jì)方法中，涉及到相應(yīng)參數(shù)的設(shè)置問題(例如LLKE中參數(shù)c的取值對估計(jì)曲線的平滑度有很大影響)。針對不同的參數(shù)設(shè)置，本文將通過仿真以及案例進(jìn)一步探討樣本量大小、估計(jì)方法對序數(shù)輪廓控制圖的IC和OC性能的影響，直觀地展現(xiàn)出選擇合適的樣本量和模型估計(jì)方法的重要性。

1 序數(shù)輪廓的監(jiān)控方法

1.1 序數(shù)輪廓的回歸模型

為了簡化表示，本文考慮的模型中僅有1個協(xié)變量。假設(shè)隨著時間的推移，收集到第j(j≥1)個輪廓。yji是第j個輪廓中的第i個響應(yīng)觀測值，xi是相應(yīng)的解釋變量(i=1,2,…，Nj)，其中Nj是第j個輪廓中的設(shè)計(jì)點(diǎn)數(shù)。假設(shè)yji是具有K個屬性級別的有序響應(yīng)變量，且yji服從參數(shù)為nji和pji的多項(xiàng)式分布，其中nji=yji1+yji2+…+yjiK，pji=(pji1,pji2,…,pjiK)T。對于不同的j和i，nji取相同的值n，yjik表示第j個輪廓中處在k級的第i個響應(yīng)變量，pjik表示yji處在k級的概率。令cjik為k級的累積概率(k=1,2，…，K)。因此cjik=Pr(yji

在對序數(shù)響應(yīng)變量和解釋變量之間的關(guān)系進(jìn)行建模時，MCCULLAGH[20]提出的比例優(yōu)勢比模型應(yīng)用較多，IZADBAKHSH等[21]將比例優(yōu)勢比模型作為擬合序數(shù)響應(yīng)變量與解釋變量之間函數(shù)關(guān)系的基礎(chǔ)，DING等[6]使用以下參數(shù)比例優(yōu)勢比模型對序數(shù)輪廓進(jìn)行建模：

logit(ck)=αk+XTβk=1,2,…,K-1

本文使用以下非參數(shù)比例優(yōu)勢比模型擬合序數(shù)輪廓的回歸模型：

logit(cjik)=gk(xjik)

(1)

j=1,2,…i=1,2,…,Njk=1,2,…,K-1

其中，gk(·)是未知的平滑函數(shù);xjik為第j個輪廓中處于k級的第i個響應(yīng)變量對應(yīng)的解釋變量。

一般來說，輪廓設(shè)計(jì)點(diǎn)可以是固定不變的(即固定設(shè)計(jì))也可以是隨機(jī)的(即隨機(jī)設(shè)計(jì))。本文只考慮固定設(shè)計(jì)，將隨機(jī)設(shè)計(jì)的方法留待以后研究。也就是說，對于不同的j，取Nj為相同的值N。對于不同的輪廓，解釋變量x=(xji,xj(i+1),…,xjN)T是已知且固定的，其中xji表示第j個輪廓的第i個響應(yīng)變量對應(yīng)的解釋變量。

序數(shù)輪廓的響應(yīng)變量yji服從多項(xiàng)式分布，根據(jù)非參數(shù)比例優(yōu)勢比模型(式(1))，可以得到第j個輪廓的對數(shù)似然函數(shù)(省略常數(shù)項(xiàng))：

(2)

式(1)的監(jiān)控問題可以看作是回歸函數(shù)gk(·)的假設(shè)檢驗(yàn)問題，原假設(shè)和備擇假設(shè)分別如下：

1.2 模型估計(jì)方法

本文主要通過3種模型估計(jì)方法對非參數(shù)比例優(yōu)勢比模型(式(1))進(jìn)行估計(jì)，分別是非參數(shù)方法LLKE、樣條法，以及參數(shù)方法Newton Raphson。對于這兩類模型估計(jì)方法，參數(shù)方法將輪廓函數(shù)擬合為線性、非線性和其他復(fù)雜模型，非參數(shù)方法則不對模型做過多的假設(shè)和約束。本文關(guān)注非參數(shù)和參數(shù)方法在模型估計(jì)中的擬合精度問題。

1.2.1局部線性核估計(jì)法

局部平滑法是估計(jì)非參數(shù)函數(shù)的常用方法之一。局部擬合方法包括局部常數(shù)擬合方法和局部多項(xiàng)式擬合方法。與前者相比，后者在擬合過程中沒有邊界效應(yīng)，擬合結(jié)果更好。本文利用局部多項(xiàng)式擬合的一種特殊情況——局部線性擬合來估計(jì)未知函數(shù)gk。

特別地，在給定的x0的鄰域內(nèi)，函數(shù)gk(xjik)通過泰勒展開式展開為線性函數(shù)：

據(jù)FAN等[22]分析可知，對數(shù)似然函數(shù)(式(2))的核函數(shù)為Kh(xji-x0)，其中h為帶寬參數(shù)。因此局部加權(quán)對數(shù)似然函數(shù)可表示為

(3)

βk的估計(jì)可通過迭代加權(quán)最小二乘法(iterative weighted least squares，IWLS)的迭代過程來獲得。YEH等[23]簡要說明了使用IWLS法估計(jì)輪廓參數(shù)的步驟。為了便于表達(dá)而不影響理解，以下公式省略了一些下標(biāo)。增廣因變量可定義為

zk=(z1k,z2k,…，zNk)T

gk=(gk(x1),gk(x2),…，gk(xN))T

Wk=diag(w1k,w2k,…，wNk)

yk=(y1k,y2k,…,yNk)T

μk=(μk(x1),μk(x2),…,μk(xN))T

此外，根據(jù)多項(xiàng)式分布的均值和方差公式以及回歸函數(shù)gk(xi)的估計(jì)，可以得到

wik=npik(1-pik)μk(xi)=npik

i=1,2，…，Nk=1,2，…,K-1

基于上述條件，可以通過求解以下IWLS方程來獲得βk的估計(jì)值：

Kh(xN-x0))

使用局部擬合方法的一個重要問題就是確定局部鄰域，該鄰域主要是由核函數(shù)Kh(·)和帶寬參數(shù)h決定的。本文選取Epanechnikov核函數(shù)：

其中，I(·)是指示函數(shù)。對于帶寬h的選擇，建議使用以下經(jīng)驗(yàn)帶寬公式[24]：

(4)

1.2.2樣條法

樣條法不再將樣本數(shù)據(jù)集當(dāng)作一個整體，而是將它劃分成一個個連續(xù)區(qū)間，劃分的點(diǎn)稱為節(jié)點(diǎn)(knot)，并用單獨(dú)的模型(線性函數(shù)或低階多項(xiàng)式函數(shù))來擬合。

樣條估計(jì)包括光滑樣條估計(jì)、多項(xiàng)式樣條估計(jì)、懲罰樣條估計(jì)和局部自適應(yīng)回歸樣條估計(jì)等方法?？紤]到本文輪廓解釋變量中N的選擇問題，以及多項(xiàng)式樣條對節(jié)點(diǎn)的個數(shù)以及位置較敏感，本文選擇多項(xiàng)式樣條對回歸函數(shù)gk進(jìn)行估計(jì)。通常情況下，當(dāng)多項(xiàng)式的階數(shù)達(dá)到3時樣條法就可獲得較高的估計(jì)精度，因此針對一個輪廓下的gk，本文使用三次多項(xiàng)式進(jìn)行擬合。

其中,θs為系數(shù);s為函數(shù)f(x)的項(xiàng)數(shù)。令gk(x)=fk(x)+ε(ε為誤差項(xiàng))，則回歸函數(shù)gk(x)的估計(jì)問題轉(zhuǎn)化為gk(x)-fk(x)=ε的最小化問題:

令

以上為使用三次樣條估計(jì)gk的大致過程，在仿真部分本文使用的是三次樣條插值法，可直接調(diào)用Python中程序包完成模型估計(jì)。

1.2.3NewtonRaphson方法

Newton Raphson方法在統(tǒng)計(jì)中廣泛應(yīng)用于求解最大似然估計(jì)。對于序數(shù)輪廓，可強(qiáng)制使用廣義線性模型擬合，然后采用Newton Raphson算法進(jìn)行參數(shù)估計(jì)[6]。本文目標(biāo)函數(shù)為對數(shù)似然函數(shù)(式(2))極大化，可轉(zhuǎn)化為求解函數(shù)l(·)的一階導(dǎo)函數(shù)l′(·)=0的問題。具體迭代步驟如下：

(1)將目標(biāo)函數(shù)l(·)在給定的x0的鄰域內(nèi)，通過泰勒展開式展開到二階形式：

(5)

當(dāng)且僅當(dāng)Δx=xji-x0無線趨近于0時，式(5)得以成立，此時式(5)等價于：

(2)計(jì)算迭代公式：

其中，v為迭代次數(shù)，v≥0。

(3)結(jié)合對數(shù)似然函數(shù)(式(2))，令

k=1,2,…,K-1

其中，Ak為關(guān)于gk的函數(shù)，則式(2)可以轉(zhuǎn)化為有關(guān)gk的K元函數(shù)f(·)：

(4)對函數(shù)f(·)求偏導(dǎo)，計(jì)算迭代公式：

Xv+1=Xv-[Hf(Xv)]-1f(Xv)

其中，H為Hessian矩陣。

(5)將Xv+1代入步驟(3)中的K元函數(shù)f(·)，重復(fù)步驟(4)、步驟(5)，直至|Xv+1-Xv|<ζ(ζ為搜索精度)，終止計(jì)算。

1.3 監(jiān)控序數(shù)輪廓的控制圖

GLR控制圖的應(yīng)用范圍較廣，例如監(jiān)控正態(tài)分布的均值或方差等。由于GLR統(tǒng)計(jì)量的良好性質(zhì)，本文基于GLR統(tǒng)計(jì)量構(gòu)建控制圖對非參數(shù)比例優(yōu)勢比模型(式(1))進(jìn)行監(jiān)控。3種估計(jì)方法下GLR統(tǒng)計(jì)量定義如下。

LLKE方法下的統(tǒng)計(jì)量：

樣條和Newton Raphson方法下的統(tǒng)計(jì)量為

控制限L(1)和L(2)可用二分法搜索，具體迭代步驟如下：

(1)在第i次迭代計(jì)算中，在區(qū)間[L(i),U(i)]中搜索L*。當(dāng)i=1時，令L(1)=0和U(1)=U，其中U是滿足以下條件的上限，即L*=U時控制圖的平均運(yùn)行鏈長(average run length，ARL)值大于預(yù)先指定的受控狀態(tài)下的ARL值A(chǔ)RL0。

(2)將搜索區(qū)間折半，即L*=l(i)=(L(i)+U(i))/2，然后計(jì)算對應(yīng)控制限L*下的運(yùn)行鏈長(run length，RL)。

2 控制圖IC性能分析

常用的階段Ⅱ控制圖性能評價指標(biāo)是ARL值?？刂茍D在IC狀態(tài)下的性能主要是通過ARL值A(chǔ)RL0進(jìn)行評估，對于OC性能，ARL值的計(jì)算存在零態(tài)和穩(wěn)態(tài)兩種假設(shè)。本文使用穩(wěn)態(tài)OC下的ARL值A(chǔ)RL1作為一個指標(biāo)來比較不同估計(jì)方法下的控制圖的OC性能。理想情況下，在ARL0相等的情況下，對于給定的偏移，具有較小ARL1的控制圖通常被認(rèn)為具有更好的性能。除此之外，用運(yùn)行鏈長的標(biāo)準(zhǔn)差(standard deviation of run length，SDRL)和運(yùn)行鏈長分布的百分位數(shù)來更全面地評估控制圖的性能。在分析控制圖的IC性能時，本文主要通過ARL0、SDRL0、運(yùn)行鏈長的10%分位數(shù)RL10%、運(yùn)行鏈長的中位數(shù)(median run length)MRL0和90%分位數(shù)RL90%來評價控制圖的性能。通過均方誤差(mean square error，MSE)來評價不同估計(jì)方法的精確度。OC狀態(tài)下的評價指標(biāo)包括ARL1、SDRL1以及MRL1指標(biāo)。

采用二分搜索法，通過5000次重復(fù)仿真，逼近相應(yīng)控制圖的控制限。為了便于說明，變化點(diǎn)設(shè)置為τ=25。為了評估每個控制圖的穩(wěn)態(tài)ARL值A(chǔ)RL1，在第(τ+1)個觀測值之前發(fā)生警報(bào)的任何仿真序列將被舍棄。對于多項(xiàng)式分布中的參數(shù)，假設(shè)n=50。另外，帶寬h由式(4)確定，常數(shù)c取3個不同的值，分別為1、1.5和2。不失一般性，假設(shè)每個序數(shù)輪廓樣本中有N=50個等距設(shè)計(jì)點(diǎn)，有

假設(shè)序數(shù)響應(yīng)變量具有3個水平，即K=3，非參數(shù)比例優(yōu)勢比模型(式(1))可以寫為

logit(cji1)=g1(xji)logit(cji2)=g2(xji)

i=1,2,…,Nj=1,2,…

本文主要考慮了廣義線性模型(generalized linear model，GLM)為受控模型的情況。

GLM受控模型為

logit(cjik)=αk+βxi

i=1,2，…,Nj=1,2,…k=1,2

如DING等[6]所述，響應(yīng)變量y的序數(shù)水平主要由截距決定，系數(shù)β反映x對y的影響。本文仍采用與DING等[6]相同的參數(shù)設(shè)置。考慮到響應(yīng)變量之間的序數(shù)信息，系數(shù)選擇為(α1,α2)=(-1,0.5)，β=0.2。

表1、表2、表3所示為在不同樣本量m下，分別使用樣條、LLKE和Newton Raphson估計(jì)方法時，GLR控制圖的ARL0、SDRL0、RL10%、MRL0和RL90%值。注意，表1、表2和表3的最后一列給出的是m=∞時控制圖的IC性能，即GLR控制圖在回歸模型已知條件下的性能。在每張表的最后一行，給出了不同樣本量下估計(jì)的均方誤差MSE。

表1 樣條方法下GLR控制圖的IC性能

一方面，從表1中可以看到隨著樣本量m的增大，GLR控制圖的各項(xiàng)評價指標(biāo)越來越接近模型已知條件下的各指標(biāo)值。在樣本量m≤800時，ARL0、SDRL0、MRL0和RL90%值呈現(xiàn)出遞增的趨勢，即樣條方法下產(chǎn)生的虛報(bào)警次數(shù)在減少。當(dāng)樣本量增至800后，再繼續(xù)增大樣本量，GLR控制圖性能的各項(xiàng)指標(biāo)不再發(fā)生變化。換言之，在樣本量m=800時，使用樣條方法就可以達(dá)到較好的估計(jì)效果，即此時樣條估計(jì)對控制圖性能的影響較小。

表2 LLKE方法下GLR控制圖的IC性能

表3 Newton Raphson 方法下GLR控制圖的IC性能

另一方面，根據(jù)MSE不難發(fā)現(xiàn)，隨著樣本量m的增大，MSE呈現(xiàn)出遞減的趨勢，即估計(jì)精度提高。最后，從表1中可以發(fā)現(xiàn)，樣本量m從30增大至800時，ARL0、SDRL0等指標(biāo)的數(shù)值跨度較大，這也說明使用樣條方法估計(jì)回歸模型時，樣本量對GLR控制圖的性能影響顯著。

表2總結(jié)了LLKE方法下GLR控制圖在不同樣本量及參數(shù)c下的IC性能。在c=2和1.5時，GLR控制圖的ARL0、SDRL0、RL10%、MRL0和RL90%指標(biāo)均隨著樣本量m的增大而遞減。當(dāng)c=2時，在樣本量m=800處GLR控制圖的IC性能與m=∞時的性能最接近，此時使用LLKE估計(jì)非參數(shù)比例優(yōu)勢比模型(式(1))對GLR控制圖的IC性能影響最小。當(dāng)m>800時，LLKE對控制圖IC性能的影響不再隨著樣本量的增大而變化。當(dāng)c=1.5時，在樣本量m=100處GLR控制圖的IC性能受LLKE估計(jì)的影響最小。在c=1的情況下，當(dāng)樣本量達(dá)到500時，LLKE估計(jì)方法下GLR控制圖的IC性能與模型已知情況下的差距最小，繼續(xù)增大樣本量，控制圖IC性能指標(biāo)不再發(fā)生變化，即m=500時，使用LLKE可獲得較好的估計(jì)效果。

觀察不同參數(shù)c下各樣本量對應(yīng)的MSE指標(biāo)，可以發(fā)現(xiàn)MSE指標(biāo)均隨著樣本量的增大而遞減，即樣本量越大，取得的估計(jì)效果越好。但是，當(dāng)樣本量較大時，MSE數(shù)值變化較小，因此，在使用LLKE方法估計(jì)回歸模型時，使用小樣本就可達(dá)到較高的估計(jì)精度。對比不同參數(shù)c下的GLR控制圖的IC指標(biāo)不難發(fā)現(xiàn)，在m≤200時，隨著c的減小，各指標(biāo)呈現(xiàn)出遞減的趨勢。當(dāng)m≥500時，各參數(shù)c下GLR控制圖的IC性能幾乎不再受樣本量大小的影響，小樣本量時控制圖已可以獲得較好的IC性能。

觀察表3，首先，隨著樣本量m的增大，在Newton Raphson 方法下GLR控制圖的各項(xiàng)IC指標(biāo)呈遞增的趨勢。當(dāng)樣本量m=200時，控制圖的性能指標(biāo)與模型已知時的性能差距較小。后續(xù)隨著樣本量繼續(xù)增大，ARL0、SDRL0、MRL0和RL90%等指標(biāo)逐漸遠(yuǎn)離m=∞時的數(shù)值，且當(dāng)m≥500時，GLR控制圖的各性能指標(biāo)不再發(fā)生變化。再者，對比不同樣本量下的MSE指標(biāo)可以發(fā)現(xiàn)，隨著樣本量的增大，MSE的遞減幅度較小，這表明在樣本量m=200后再繼續(xù)增大樣本數(shù)量對Newton Raphson方法的估計(jì)效果影響不大。因此，在使用Newton Raphson方法估計(jì)GLM時，小樣本量也可達(dá)到較好的估計(jì)效果。

綜合表1、表2和表3，對比3種模型估計(jì)方法所需的最小最優(yōu)樣本量，在IC模型是GLM時，與非參數(shù)方法相比，使用參數(shù)方法Newton Raphson可以達(dá)到較高的估計(jì)精度，對階段Ⅱ的GLR控制圖的IC性能影響較小。

由以上分析可知，模型估計(jì)會影響GLR控制圖的IC性能，而適當(dāng)增大樣本量會降低該影響，但由于成本或其他限制，收集大量樣本比較困難，在這種情況下，調(diào)整控制限也可以補(bǔ)償模型估計(jì)帶來的影響。本文通過模擬實(shí)驗(yàn)，調(diào)整GLR控制圖在不同估計(jì)方法下的控制限，使其具有預(yù)定的ARL0約等于200。樣條方法下對應(yīng)的GLR控制圖的調(diào)整后的控制限見表4。根據(jù)表4可知，在樣本量一定的情況下，控制圖各性能指標(biāo)相對調(diào)整控制限之前更加接近m=∞時的性能指標(biāo)。當(dāng)m=30時，僅看MRL0指標(biāo)，表1中調(diào)整控制限之前MRL0為38，表4中MRL0為124，顯然后者更接近m=∞時的MRL0值122。在LLKE和Newton Raphson方法下，可以得到相同的結(jié)論。

表4 調(diào)整控制限后樣條方法下GLR控制圖的IC性能

3 控制圖OC性能分析

對于OC模型，本文僅選取了一種有代表性且易于理解的OC模型進(jìn)行仿真研究：只改變模型系數(shù)，響應(yīng)變量與解釋變量之間的函數(shù)關(guān)系結(jié)構(gòu)不變。失控模型如下：

logit(cjik)=(αk+δ1)+(β+δ2)xi

i=1,2,…,Nj=1,2,…k=1,2

為了評估失控狀態(tài)下使用不同的估計(jì)方法對階段Ⅱ的GLR控制圖OC性能的影響，本文主要考慮了3個樣本量范圍。為了公平比較控制圖的OC性能，應(yīng)確保ARL0相同，因此采用調(diào)整后的控制限。首先在樣本量一定的情況下，評價樣條和Newton Raphson方法對GLR控制圖OC性能的影響。接著，當(dāng)樣本量發(fā)生變化時，分析控制圖性能指標(biāo)的變化情況。最后，在不同樣本量下，比較樣條、LLKE和Newton Raphson方法下控制圖的OC性能指標(biāo)，找到在IC模型為GLM時，參數(shù)和非參數(shù)模型估計(jì)方法適用的樣本量范圍。

表5和表6分別總結(jié)了在使用樣條和Newton Raphson估計(jì)方法的情況下，當(dāng)參數(shù)β發(fā)生偏移時，GLR控制圖在樣本量m=70和∞下的ARL1、SDRL1和MRL1的數(shù)值。從表5和表6中可以發(fā)現(xiàn)，GLR控制圖的OC性能受到模型估計(jì)的影響。當(dāng)使用樣條估計(jì)模型時，如果參數(shù)發(fā)生的偏移較小，則該估計(jì)方法會帶來較大的負(fù)面影響；當(dāng)偏移較大時，則估計(jì)對控制圖性能的影響較小。使用Newton Raphson估計(jì)回歸模型時，無論參數(shù)β發(fā)生的偏移大小，該方法均會對GLR控制圖的OC性能產(chǎn)生明顯的負(fù)面影響。

表5 樣條估計(jì)方法對GLR控制圖OC性能的影響(β偏移)

表6 Newton Raphson估計(jì)方法對GLR控制圖OC性能的影響(β發(fā)生偏移)

圖1和圖2分別給出了樣條和Newton Raphson方法在3種不同的樣本量下，當(dāng)參數(shù)α發(fā)生偏移時，GLR控制圖的OC性能指標(biāo)。從圖1中可以看到，在同一樣本量m下，隨著偏移的增大，ARL1、SDRL1和MRL1均呈現(xiàn)出遞減的趨勢，即參數(shù)發(fā)生的偏移較大時，GLR控制圖能夠快速檢測出異常并發(fā)出警報(bào)。另外，在偏移量相同的情況下，m=30時的控制圖的各項(xiàng)OC指標(biāo)遠(yuǎn)遠(yuǎn)超過m=200時的指標(biāo)，即GLR控制圖的檢出力度隨著樣本量的增大而增強(qiáng)。因此，適當(dāng)增大樣本量可以減小模型估計(jì)對控制圖OC性能的影響。由圖2可以得出類似的結(jié)論，但與樣條方法不同的是，使用Newton Raphson估計(jì)回歸模型時，隨著樣本量m的增大，同一偏移下的各OC指標(biāo)的差距較小。因此，使用Newton Raphson方法估計(jì)模型時增大樣本量m，對控制圖檢出力的影響并不顯著。

(c)MRL1的對比圖1 樣條方法下樣本量對GLR控制圖OC性能的影響(α偏移)Fig.1 The effects of sample sizes on OC performanceof GLR control chart under spline method(α shifts)

(a)ARL1的對比 (b)SDRL1的對比

(c)MRL1的對比圖2 Newton Raphson方法下樣本量對GLR控制圖OC性能的影響(α偏移)Fig.2 The effects of sample sizes on OC performanceof GLR control chart under Newton Raphsonmethod(α shifts)

通過以上對不同樣本量下使用樣條和Newton Raphson方法時的GLR控制圖的OC性能分析可知，樣本量對階段Ⅱ控制圖性能的影響因估計(jì)方法的不同而有所不同，所以在實(shí)際應(yīng)用時，應(yīng)根據(jù)樣本量的大小謹(jǐn)慎選擇模型估計(jì)方法。在LLKE方法下，式(4)中的參數(shù)c對估計(jì)曲線的平滑度有較大影響，進(jìn)而會影響最終的模型估計(jì)結(jié)果。因此，對于LLKE方法，本文重點(diǎn)關(guān)注參數(shù)c對GLR控制圖OC性能的影響。表7所示為樣本量固定為70的情況下，使用LLKE方法估計(jì)模型時，GLR控制圖的OC性能指標(biāo)。根據(jù)表7可得，在同一參數(shù)c下，當(dāng)參數(shù)β發(fā)生偏移時，隨著偏移量的增大，GLR控制圖的OC性能指標(biāo)呈遞減的趨勢；當(dāng)偏移量固定時，隨著參數(shù)c的增大，ARL1、SDRL1和MRL1大致呈遞減的趨勢，GLR控制圖的OC性能隨著c的遞增而增強(qiáng)。

表7 LLKE估計(jì)方法對GLR控制圖OC性能的影響(β偏移)

對比分析不同的樣本量下采用這3種估計(jì)方法的控制圖OC性能指標(biāo)，得出3種模型估計(jì)方法分別適用的樣本量范圍。本文選取了樣本量分別為30,70,200的3種情況進(jìn)行簡單解釋。當(dāng)參數(shù)α偏移時，GLR控制圖的OC性能指標(biāo)數(shù)據(jù)見表8。

表8 不同樣本量下樣條、LLKE(c=1.5)和Newton Raphson三種估計(jì)方法對GLR控制圖OC性能的影響(α偏移)

由表8可以發(fā)現(xiàn)，首先，對于同一種估計(jì)方法，隨著樣本量的增大，GLR控制圖的ARL1、SDRL1和MRL1指標(biāo)均呈現(xiàn)出遞減的趨勢，即在大樣本量情況下控制圖的OC性能更優(yōu)，該結(jié)論與圖1和圖2一致。其次，當(dāng)樣本量固定在30時，比較3種估計(jì)方法下GLR控制圖的OC性能可以發(fā)現(xiàn)，Newton Raphson方法下控制圖的OC性能最優(yōu)，無論偏移大小，控制圖均能夠及時發(fā)現(xiàn)變點(diǎn)，LLKE(c=1.5)次之，而樣條方法下控制圖的OC性能最差。在m=70和200時可以得到同樣的結(jié)論。因此，當(dāng)IC模型為GLM且樣本量m≤200時，使用參數(shù)方法Newton Raphson可以獲得較好的GLR控制圖OC性能。

在仿真過程中，本文假設(shè)IC模型為GLM，該假設(shè)與參數(shù)方法Newton Raphson的假設(shè)一致，因此使用Newton Raphson估計(jì)模型時，GLR控制圖的性能較好。但對于復(fù)雜的函數(shù)關(guān)系，無法使用GLM擬合模型。通過以上仿真可以發(fā)現(xiàn)，在IC模型為GLM時，LLKE方法下的控制圖性能僅次于Newton Raphson，因此在模型為非GLM或未知時，推薦使用LLKE。

4 案例研究

為了進(jìn)一步檢驗(yàn)LLKE、樣條和Newton Rapson方法的估計(jì)效果，本文將3種方法應(yīng)用于一個汽車制造商的保修索賠的數(shù)據(jù)集。汽車是最常見的帶有保修合同的產(chǎn)品之一，在保修期內(nèi)，制造商有義務(wù)免費(fèi)維修或更換保修合同范圍內(nèi)的瑕疵產(chǎn)品。汽車投入使用后，將每天(每周或每月)發(fā)生的保修索賠次數(shù)記錄在保修數(shù)據(jù)集中，通過對保修索賠數(shù)據(jù)建模分析，可盡早發(fā)現(xiàn)質(zhì)量或現(xiàn)場可靠性問題[25]。通常假定同一天生產(chǎn)的產(chǎn)品具有相同的可靠性[26]。

在本文所使用的保修索賠數(shù)據(jù)集中，每周生產(chǎn)的汽車總數(shù)量中的累計(jì)保修索賠數(shù)量是按照汽車使用月份記錄的。為了排除其他因素的影響，本文只考慮在特定車型和型號下汽車發(fā)動機(jī)的保修索賠情況。因?yàn)榘l(fā)動機(jī)的保修期通常不少于兩年，所以選擇對汽車銷售后24個月的索賠數(shù)據(jù)進(jìn)行分析。數(shù)據(jù)集中提供了故障類型(從工程角度)和每種故障類型下的保修索賠數(shù)量。通過整理數(shù)據(jù)發(fā)現(xiàn)，所有故障類型按照發(fā)生頻次可歸納為輕微、中等、嚴(yán)重3種故障類型。案例中的質(zhì)量特征是一個輪廓，其中響應(yīng)變量是不同故障類型的累計(jì)保修索賠數(shù)量，解釋變量是使用月份。數(shù)據(jù)集中共有52個保修索賠的輪廓，圖3顯示了6個具有代表性的輪廓(第19個輪廓到第24個輪廓)的logit(c)圖。由圖3可以看出，圖3c與圖3d之間存在明顯的差異，即從第22個輪廓開始，函數(shù)模型開始發(fā)生變化。因此，在案例研究中前21個輪廓作為IC參考數(shù)據(jù)。

(a)第19個輪廓 (b)第20個輪廓

(c)第21個輪廓 (d)第22個輪廓

(e)第23個輪廓 (f)第24個輪廓圖3 IC和OC狀態(tài)下的輪廓的logit(c)曲線Fig.3 The logit(c) curve of the profile in IC andOC states

采用非參數(shù)比例優(yōu)勢比模型(式(1))來表示響應(yīng)變量y與解釋變量X之間的函數(shù)關(guān)系。

綜上，當(dāng)IC模型是GLM時，本文推薦使用Newton Raphson，該方法所需估計(jì)樣本量較少，且可以獲得較好的估計(jì)效果。當(dāng)IC模型形式不明確時，本文推薦使用非參數(shù)方法LLKE和樣條。

(a)基準(zhǔn) (b)LLKE

(c)樣條(d)Newton Raphson圖4 使用LLKE、樣條和Newton Raphson方法的模型估計(jì)結(jié)果與基準(zhǔn)對比Fig.4 Comparisons of Benchmark with the estimationresults using LLKE, spline and Newton Raphson method

5 結(jié)論及展望

本文評估了使用樣條、LLKE和Newton Raphson 3種估計(jì)方法對序數(shù)輪廓控制圖的IC和OC性能的影響。首先，本文通過ARL0、SDRL0、RL10%、MRL0和RL90%等指標(biāo)評價了3種模型估計(jì)方法對GLR控制圖IC性能的影響。通過對比發(fā)現(xiàn)，一方面增加樣本量可以減少模型估計(jì)對控制圖性能的影響，但是當(dāng)樣本量增加到一定程度之后，繼續(xù)增加樣本量所到來的效果并不理想。另一方面在IC模型為廣義線性模型(GLM)時，3種估計(jì)方法中，Newton Raphson方法對樣本量大小的要求不高。其次，當(dāng)IC模型由GLM變?yōu)榱硪环N線性模型時，本文評估了GLR控制圖在3種估計(jì)方法下的OC性能。通過ARL1、SDRL1和MRL1指標(biāo)可以發(fā)現(xiàn)，在樣本量固定的情況下，樣條、LLKE和Newton Raphson均會影響控制圖OC性能，但參數(shù)方法Newton Raphson可以獲得較好的估計(jì)效果，即在OC模型仍為GLM時Newton Raphson對GLR控制圖的OC性能影響較小。再者，本文對比了不同樣本量下使用3種估計(jì)方法時控制圖的OC指標(biāo)，發(fā)現(xiàn)當(dāng)IC模型為GLM時，參數(shù)方法較非參數(shù)方法可以達(dá)到較高的估計(jì)精度，且對控制圖性能影響較小。最后，本文通過一個實(shí)際案例進(jìn)一步說明了當(dāng)IC模型不明確時，非參數(shù)方法可以獲得較高的估計(jì)結(jié)果精確度。

本文存在以下兩點(diǎn)不足。第一，本文主要針對IC模型為GLM的情形研究了模型估計(jì)對序數(shù)輪廓控制圖性能的影響，對其他更復(fù)雜的IC模型沒有做細(xì)致的說明。在實(shí)際應(yīng)用場景中，IC模型的具體表現(xiàn)形式不得而知，因此仿真情景中，可適當(dāng)增加較復(fù)雜的模型以進(jìn)行更完善的研究。第二，LLKE中的參數(shù)c對估計(jì)精度有較大的影響[19]，關(guān)于模型估計(jì)對序數(shù)輪廓控制圖性能的影響以及LLKE方法下的參數(shù)c的影響可做進(jìn)一步研究。