曹程明, 馬義中

(南京理工大學經(jīng)濟管理學院,江蘇南京210094)

1 引 言

在高質量過程監(jiān)控中,監(jiān)控事件間隔時間是一個重要問題[1].當事件發(fā)生率不為定值時,用威布爾分布描述事件間隔時間更為合適.Woodall 等[2]、Ali 等[3]建議研究威布爾分布的TBE 控制圖.威布爾更新過程由一系列獨立同分布的威布爾分布的事件間隔時間構成,其在產(chǎn)品、設備或者生產(chǎn)系統(tǒng)的可靠性建模及監(jiān)控中廣泛應用[4].

威布爾分布的形狀參數(shù)與產(chǎn)品的失效機制緊密相關,形狀參數(shù)的變化往往意味著所監(jiān)控的產(chǎn)品或系統(tǒng)的根本性的變化,因此相比較于監(jiān)控樣本均值,監(jiān)控威布爾分布的形狀參數(shù)更有意義[4].Pascual[5]將威布爾分布轉換成最小極值分布,然后基于樣本極差得到1/β的無偏估計量作為監(jiān)控統(tǒng)計量,設計了EWMA 控制圖監(jiān)控威布爾分布的形狀參數(shù).Pascual 等[6]和Pascual 等[7]采用了相似的方法,將樣本極差作為監(jiān)控統(tǒng)計量對威布爾分布的形狀參數(shù)進行監(jiān)控.但是應用這些控制圖時,為了得到控制限需要進行復雜的計算.為了避免復雜的計算,Zhang 等[4]利用形狀參數(shù)的樞軸量構建了一個近似服從正態(tài)分布的統(tǒng)計量,并在此基礎上設計了Beta 控制圖,仿真分析結果說明Beta 控制圖能夠有效的探測到威布爾分布形狀參數(shù)的偏移.但是其研究中只應用20 個過程穩(wěn)定時的樣本數(shù)據(jù)估計過程參數(shù),且沒有分析參數(shù)估計對其所設計的控制圖性能的影響.由于參數(shù)估計量固有的波動,在同一設計方法下設計出的控制圖的性能將呈現(xiàn)差異,這稱之為使用者與使用者之間的波動.這種波動使得使用者不能確定所設計控制圖的錯誤警報概率,影響其對控制圖的信心.

在應用控制圖的第一階段,收集大量的樣本數(shù)據(jù)通常是不現(xiàn)實的,貝葉斯估計利用了過程的先驗信息,可以減少參數(shù)估計對控制圖性能的影響[8].Raubenheimer 等[9]設計了Bayesian c 圖,分析結果表明Bayesian c 圖具有良好的性能.Kumar 等[10]設計了Bayesian-Tr 控制圖,分析結果表明Bayesian-Tr 控制圖的性能優(yōu)于傳統(tǒng)的Tr 控制圖.Lee 等[11]設計了基于貝葉斯理論的Bernoulli CUSUM 控制圖,提高了控制圖的性能.

針對威布爾更新過程,考慮威布爾分布形狀參數(shù)偏移的探測問題以及實際中樣本數(shù)據(jù)較少的問題,本文在Beta 控制圖的基礎上,結合貝葉斯理論,設計了Bayesian-Beta 控制圖監(jiān)控威布爾更新過程的形狀參數(shù).首先,給出了Bayesian-Beta 控制圖的統(tǒng)計設計方法及其監(jiān)控程序;其次,給出了過程參數(shù)的貝葉斯估計方法;然后,應用AARL 和SDARL 及蒙特卡洛仿真分析了參數(shù)估計對控制圖性能的影響,結果表明Bayesian-Beta控制圖受參數(shù)估計量波動的影響更小,所提出的貝葉斯估計方法中先驗信息的誤差對控制圖性能的影響小,可操作性強;最后,用實際案例說明Bayesian-Beta 控制圖在實際中的應用.

2 Bayesian-Beta 控制圖的設計

2.1 監(jiān)控統(tǒng)計量

兩參數(shù)威布爾分布的密度函數(shù)如下

其中α,β >0,分別是威布爾分布的尺度和形狀參數(shù).

令{x1,x2,...,xn}為威布爾分布的一組樣本值,則α和β的極大似然估計α和β可由下式得到,即

Zhang 等[4]的仿真實驗表明,λ=?0.5 時,近似服從正態(tài)分布.的均值Ec[]和方差Varc()分別為其中En=E[(/β)λ],Vn=Var[(/β)λ].(/β)λ是樞軸量,因此En和Vn可通過蒙特卡洛仿真得到.En和Vn的參考值可參照文獻[4].

2.2 控制限的設計

對于一系列來自威布爾更新過程的事件間隔時間的數(shù)據(jù),設定窗口大小為n,{x1,x2,...,xn}為第一個樣本,{x2,x3,...,xn+1}為第二個樣本,依此類推,根據(jù)式(3)可計算統(tǒng)計量S的值.由2.1 節(jié)可知統(tǒng)計量服從均值為的正態(tài)分布.Zhang 等[4]設計了Beta 控制圖,其上控制限UCL和中心線CL 以及下控制限LCL 分別為

其中n為計算監(jiān)控統(tǒng)計量的窗口大小,其中k為控制圖系數(shù).

顯然控制圖的ARL0值由n和k決定,不同n和k的ARL0值可參照Zhang 等[4].βic為過程受控時的β值.在應用控制圖的第一階段,收集大量的樣本數(shù)據(jù)通常是不現(xiàn)實的,為了利用過程的先驗信息,本文使用貝葉斯方法估計β值,設計Bayesian-Beta 控制圖,其控制限的設計與Beta 控制圖相同.

2.3 控制圖的監(jiān)控程序

質量控制的第一階段: 確定過程的受控狀態(tài).首先根據(jù)第一階段的m個樣本數(shù)據(jù)得到參數(shù)β的貝葉斯估計值;然后選取合適的窗口大小n和系數(shù)k計算控制限、打點,判斷過程是否受控,若有異常點且異常點是由異常原因造成,則需要查找異常原因并消除,調(diào)整生產(chǎn)參數(shù),采集數(shù)據(jù)重新估計參數(shù)β,計算控制限和打點,直到確定過程處于受控狀態(tài),將此時的作為過程參數(shù)的估計值,開始第二階段的監(jiān)控.第二階段:監(jiān)控階段.計算,打點,若打點值超出控制限,則可判斷過程失控,查找失控原因,采取措施使過程回到受控狀態(tài)(當i

3 貝葉斯估計

3.1 參數(shù)估計的步驟

式(4)是建立在已知參數(shù)βic的基礎上，在實際中參數(shù)βic通常是未知的.假設X=(x1,x2,...,xm)為第一階段的樣本數(shù)據(jù).

令θ=1/αβ,則式(1)變?yōu)?/p>

假設參數(shù)θ的先驗分布為伽馬分布Gamma(c,d),c,d>0,β的先驗分布為伽馬分布Gamma(e,f),e,f >0 且獨立于θ[13].可得θ和β的聯(lián)合后驗密度為

因此,θ和β的后驗密度為

在平方損失下,β的點估計為其后驗分布的均值.這里采用獨立鏈M-H 算法從其后驗密度中對β進行采樣,具體步驟如下:

步驟1任意選擇一個初始值β(i),并令i=0;

步驟2從伽馬分布Gamma(e,f)中生成一個備選值β?,并從均勻分布U(0,1)中抽取u;

步驟3若u≤a(β(i),β?),令β(i+1)=β?,否則令β(i+1)=βi,接受概率為

步驟4令i ←i+1,并返回到步驟2,重復步驟2、步驟3,M=10 000 次.

則參數(shù)β的估計值為

3.2 先驗分布參數(shù)的取值

令prior(β)和prior(α)分別為β和α的先驗值,MLE(β)和MLE(α)分別為參數(shù)β和α在樣本X=(x1,x2,...,xm)下的極大似然估計.則β和θ先驗分布的參數(shù)值可由式(5)和式(6)得到

4 控制圖的性能分析

4.1 受控狀態(tài)下控制圖的性能分析

過程受控時,為了分析參數(shù)估計對控制圖性能的影響,考慮當?shù)谝浑A段的樣本數(shù)m為20,50,100,200,250 和500, 參數(shù)β真實值分別為0.5, 1 和2 時控制圖的性能.不失一般性, 令α= 1.表1、表2 和表3分別是n= 11,k= 2.2;n= 15,k= 2.1 和n= 20,k= 2.0 即ARL0= 100 時Bayesian-Beta 控制圖的AARL 和SDARL 值.為了比較Bayesian-Beta 控制圖(貝葉斯估計)和Beta 控制圖(極大似然估計)的性能差異,表1、表2 和表3 同時給出了Beta 控制圖的AARL 和SDARL 值.具體的仿真步驟如下:

步驟1初始化m,n,k,β,α,并令先驗值為prior(β)=β和prior(α)=α;

步驟2生成m個第一階段的樣本數(shù)據(jù),估計過程參數(shù)后進入過程監(jiān)控階段;

步驟3產(chǎn)生一系列監(jiān)控階段的數(shù)據(jù),計算并打點,當點超出控制限時,認為過程失控,此時的i為運行鏈長,;

步驟4重復步驟3 共10 000 次,計算i的均值得到ARL;

步驟5重復步驟2 到步驟4 共2 000次,計算ARL 的均值和標準差得AARL 和SDARL.

從表1、表2和表3 中可以看出隨著樣本數(shù)m的增大,Bayesian-Beta 控制圖和Beta 控制圖的AARL 逐漸接近100,SDARL 越來越小,即隨著樣本數(shù)m的增大,參數(shù)未知的Bayesian-Beta 控制圖和Beta 控制圖的性能越來越接近參數(shù)已知的情況.對比Bayesian-Beta 控制圖和Beta 控制圖的性能可以發(fā)現(xiàn)Bayesian-Beta控制圖的AARL 更大,SDARL 更小,例如,當m=20,β=0.5 時,Bayesian-Beta 控制圖的AARL=75.73,SDARL = 26.53, Beta 控制圖的AARL = 63.63, SDARL = 30.79, 說明用貝葉斯估計可以有效降低參數(shù)估計對控制圖性能的影響.Zhang 等[14]建議SDARL 值應小于設定值ARL0的10%,對于Bayesian-Beta控制圖, 當樣本數(shù)m為200 時, 其SDARL 值在設定值ARL0的10%左右, 而Beta 控制圖, 當樣本數(shù)m為250 甚至500 時,其SDARL 值才在設定值ARL0的10%左右.

對比表1、表2 和表3 可以看出, 當ARL0相同, 窗口大小n不同時控制圖的性能略有差異.當樣本數(shù)和參數(shù)β相同如m= 20,β= 0.5 時,n= 11 時的AARL = 75.73, SDARL = 25.63,n= 15 時的AARL = 71.74,SDARL = 26.57,n= 20 時的AARL = 66.90,SDARL = 27.15,說明參數(shù)未知情況下,窗口大小n較小時受控狀態(tài)下的控制圖的性能更穩(wěn)定.

表1 受控狀態(tài)下控制圖的性能(n=11,k =2.2)Table 1 Control chart performance in controlled state(n=11,k =2.2)

表2 受控狀態(tài)下控制圖的性能(n=15,k =2.1)Table 2 Control chart performance in controlled state(n=15,k =2.1)

表3 受控狀態(tài)下控制圖的性能(n=20,k =2.0)Table 3 Control chart performance in controlled state(n=20,k =2.0)

4.2 先驗信息誤差靈敏度的分析

在貝葉斯估計中, 先驗信息對估計結果有一定的影響.在確定參數(shù)β和θ先驗分布的參數(shù)值時利用了參數(shù)β和α的先驗值prior(β)和prior(α), 顯然prior(β)和prior(α)會影響到控制圖的性能.令δ1= prior(β)/β,δ2= prior(α)/α分別為參數(shù)β和α的先驗值prior(β)和prior(α)和其真實值的誤差系數(shù).因為δ1或δ2的增大(或減小)對控制圖性能的影響相同, 所以分析中考慮δ1或δ2取相同值的情況.表4 給出了誤差系數(shù)δ1或δ2為(1.15,1.10,1.05,1.00,0.95,0.90,0.85)時的控制圖的性能,其中n= 11,k=2.2,β=1,α=1.

表4 不同誤差系數(shù)下控制圖的性能對比Table 4 Performance comparison of control charts with different error coefficients

當誤差系數(shù)δ1=δ2>1 時, 即參數(shù)β和α的先驗值prior(β)和prior(α)大于其真實值時, 控制圖的性能隨著誤差系數(shù)的增大越來越差, 但仍然比Beta 控制圖的性能更好, 例如m= 20 時, 誤差系數(shù)δ1=δ2= 1.15 的AARL = 67.51, SDARL = 25.66,δ1=δ2= 1.10 的AARL = 70.82,SDARL = 24.60,δ1=δ2= 1.05 的AARL = 73.58, SDARL = 23.32, Beta 控制圖的AARL = 61.47,SDARL = 29.14.當誤差系數(shù)δ1=δ2<1 時,即參數(shù)β和α的先驗值prior(β)和prior(α)小于其真實值時,控制圖的性能隨著誤差系數(shù)的減小越來越好.

綜上, 采用貝葉斯方法估計參數(shù)能夠有效減少參數(shù)估計量的波動對控制圖性能的影響, 即使用者與使用者之間的波動問題, 增強了使用者對控制圖的信心; 先驗信息值的誤差在真實值的15% 以內(nèi),Bayesian-Beta 控制圖都有著較好的性能,說明所提貝葉斯估計方法具有可操作性.

5 案例分析

表5 的數(shù)據(jù)是某生產(chǎn)系統(tǒng)中生產(chǎn)出來的硬盤驅動器的加速壽命試驗的失效數(shù)據(jù)[4].取前20 個數(shù)據(jù)為第一階段的樣本數(shù)據(jù),n= 11,k= 2.2 即ARL0= 100,計算得βic= 1.238,LCL = 0.624,CL = 0.867,UCL=1.110.剩余的數(shù)據(jù)作為第二階段的監(jiān)控數(shù)據(jù),監(jiān)控圖如圖1 所示.

表5 硬盤驅動器失效時間Table 5 failure times of hard disk drives

圖1 中的實心圓點表示第一階段的打點值,可以看出第一階段過程是穩(wěn)定的;其余點是第二階段的打點值,超出控制限的打點值用方框實心點表示.圖中有一個明顯的上升的趨勢,說明形狀參數(shù)的值在減小,意味著該硬盤驅動的可靠性在提高.事實上該硬盤驅動器的可靠性確實提高了,圖中的上升的趨勢符合事實.可見Bayesian-Beta 控制圖對過程偏移具有較好的探測能力.

圖1 某硬盤驅動器的生產(chǎn)過程監(jiān)控圖Fig.1 The control chart of a hard disk drives’production process

6 結束語

針對威布爾更新過程,考慮威布爾分布形狀參數(shù)偏移的探測問題以及實際中樣本數(shù)據(jù)較少的問題,設計了Bayesian-Beta 控制圖,給出了控制圖的監(jiān)控程序和參數(shù)的貝葉斯估計方法.相比較于Beta 控制圖,在過程受控時Bayesian-Beta 控制圖的AARL 更接近于ARL0,SDARL 更小,即Bayesian-Beta 控制圖的性能受參數(shù)估計量波動的影響更小.所提出的貝葉斯估計方法中先驗信息的誤差對控制圖性能的影響小,可操作性強.案例分析表明Bayesian-Beta 控制圖在實際應用中能夠有效探測到過程的偏移.

由于監(jiān)控統(tǒng)計量之間的相關性,本文未能給出ARL 的理論計算方法,另外如何處理數(shù)據(jù)之間的相關性是一個需要進一步研究的問題.