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

曹程明, 馬義中

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

1 引 言

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

威布爾分布的形狀參數(shù)與產(chǎn)品的失效機制緊密相關(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ù)估計量波動的影響更小,所提出的貝葉斯估計方法中先驗信息的誤差對控制圖性能的影響小,可操作性強;最后,用實際案例說明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ào)整生產(chǎn)參數(shù),采集數(shù)據(jù)重新估計參數(shù)β,計算控制限和打點,直到確定過程處于受控狀態(tài),將此時的作為過程參數(shù)的估計值,開始第二階段的監(jiān)控.第二階段:監(jiā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 且獨立于θ[13].可得θ和β的聯(lián)合后驗密度為

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

在平方損失下,β的點估計為其后驗分布的均值.這里采用獨立鏈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ù),計算并打點,當(dāng)點超出控制限時,認(rèn)為過程失控,此時的i為運行鏈長,;

步驟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ù)估計量的波動對控制圖性能的影響, 即使用者與使用者之間的波動問題, 增強了使用者對控制圖的信心; 先驗信息值的誤差在真實值的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 中的實心圓點表示第一階段的打點值,可以看出第一階段過程是穩(wěn)定的;其余點是第二階段的打點值,超出控制限的打點值用方框?qū)嵭狞c表示.圖中有一個明顯的上升的趨勢,說明形狀參數(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ù)估計量波動的影響更小.所提出的貝葉斯估計方法中先驗信息的誤差對控制圖性能的影響小,可操作性強.案例分析表明Bayesian-Beta 控制圖在實際應(yīng)用中能夠有效探測到過程的偏移.

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