基于貝葉斯估計的威布爾更新過程統(tǒng)計控制

曹程明, 馬義中

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

1 引 言

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

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

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

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

2 Bayesian-Beta 控制圖的設(shè)計

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的參考值可參照文獻(xiàn)[4].

2.2 控制限的設(shè)計

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

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

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

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

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

3 貝葉斯估計

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

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

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

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

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

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

步驟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,重復(fù)步驟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ù)估計對控制圖性能的影響,考慮當(dāng)?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ù)后進(jìn)入過程監(jiān)控階段;

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

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

步驟5重復(fù)步驟2 到步驟4 共2 000次,計算ARL 的均值和標(biāo)準(zhǔn)差得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 更小,例如,當(dāng)m=20,β=0.5 時,Bayesian-Beta 控制圖的AARL=75.73,SDARL = 26.53, Beta 控制圖的AARL = 63.63, SDARL = 30.79, 說明用貝葉斯估計可以有效降低參數(shù)估計對控制圖性能的影響.Zhang 等[14]建議SDARL 值應(yīng)小于設(shè)定值A(chǔ)RL0的10%,對于Bayesian-Beta控制圖, 當(dāng)樣本數(shù)m為200 時, 其SDARL 值在設(shè)定值A(chǔ)RL0的10%左右, 而Beta 控制圖, 當(dāng)樣本數(shù)m為250 甚至500 時,其SDARL 值才在設(shè)定值A(chǔ)RL0的10%左右.

對比表1、表2 和表3 可以看出, 當(dāng)ARL0相同, 窗口大小n不同時控制圖的性能略有差異.當(dāng)樣本數(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 先驗信息誤差靈敏度的分析

在貝葉斯估計中, 先驗信息對估計結(jié)果有一定的影響.在確定參數(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

當(dāng)誤差系數(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.當(dāng)誤差系數(shù)δ1=δ2<1 時,即參數(shù)β和α的先驗值prior(β)和prior(α)小于其真實值時,控制圖的性能隨著誤差系數(shù)的減小越來越好.

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

5 案例分析

表5 的數(shù)據(jù)是某生產(chǎn)系統(tǒng)中生產(chǎn)出來的硬盤驅(qū)動器的加速壽命試驗的失效數(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 硬盤驅(qū)動器失效時間Table 5 failure times of hard disk drives

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

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

6 結(jié)束語

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

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