曹程明, 馬義中

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

1 引 言

在高質(zhì)量過(guò)程監(jiān)控中,監(jiān)控事件間隔時(shí)間是一個(gè)重要問(wèn)題[1].當(dāng)事件發(fā)生率不為定值時(shí),用威布爾分布描述事件間隔時(shí)間更為合適.Woodall 等[2]、Ali 等[3]建議研究威布爾分布的TBE 控制圖.威布爾更新過(guò)程由一系列獨(dú)立同分布的威布爾分布的事件間隔時(shí)間構(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/β的無(wú)偏估計(jì)量作為監(jiān)控統(tǒng)計(jì)量,設(shè)計(jì)了EWMA 控制圖監(jiān)控威布爾分布的形狀參數(shù).Pascual 等[6]和Pascual 等[7]采用了相似的方法,將樣本極差作為監(jiān)控統(tǒng)計(jì)量對(duì)威布爾分布的形狀參數(shù)進(jìn)行監(jiān)控.但是應(yīng)用這些控制圖時(shí),為了得到控制限需要進(jìn)行復(fù)雜的計(jì)算.為了避免復(fù)雜的計(jì)算,Zhang 等[4]利用形狀參數(shù)的樞軸量構(gòu)建了一個(gè)近似服從正態(tài)分布的統(tǒng)計(jì)量,并在此基礎(chǔ)上設(shè)計(jì)了Beta 控制圖,仿真分析結(jié)果說(shuō)明Beta 控制圖能夠有效的探測(cè)到威布爾分布形狀參數(shù)的偏移.但是其研究中只應(yīng)用20 個(gè)過(guò)程穩(wěn)定時(shí)的樣本數(shù)據(jù)估計(jì)過(guò)程參數(shù),且沒(méi)有分析參數(shù)估計(jì)對(duì)其所設(shè)計(jì)的控制圖性能的影響.由于參數(shù)估計(jì)量固有的波動(dòng),在同一設(shè)計(jì)方法下設(shè)計(jì)出的控制圖的性能將呈現(xiàn)差異,這稱之為使用者與使用者之間的波動(dòng).這種波動(dòng)使得使用者不能確定所設(shè)計(jì)控制圖的錯(cuò)誤警報(bào)概率,影響其對(duì)控制圖的信心.

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

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

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

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

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

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

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

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

2.2 控制限的設(shè)計(jì)

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

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

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

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

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

3 貝葉斯估計(jì)

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

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

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

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

因此,θ和β的后驗(yàn)密度為

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

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

步驟2從伽馬分布Gamma(e,f)中生成一個(gè)備選值β?,并從均勻分布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ù)β的估計(jì)值為

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

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

4 控制圖的性能分析

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

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

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

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

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

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

步驟5重復(fù)步驟2 到步驟4 共2 000次,計(jì)算ARL 的均值和標(biāo)準(zhǔn)差得AARL 和SDARL.

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

對(duì)比表1、表2 和表3 可以看出, 當(dāng)ARL0相同, 窗口大小n不同時(shí)控制圖的性能略有差異.當(dāng)樣本數(shù)和參數(shù)β相同如m= 20,β= 0.5 時(shí),n= 11 時(shí)的AARL = 75.73, SDARL = 25.63,n= 15 時(shí)的AARL = 71.74,SDARL = 26.57,n= 20 時(shí)的AARL = 66.90,SDARL = 27.15,說(shuō)明參數(shù)未知情況下,窗口大小n較小時(shí)受控狀態(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 先驗(yàn)信息誤差靈敏度的分析

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

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

當(dāng)誤差系數(shù)δ1=δ2>1 時(shí), 即參數(shù)β和α的先驗(yàn)值prior(β)和prior(α)大于其真實(shí)值時(shí), 控制圖的性能隨著誤差系數(shù)的增大越來(lái)越差, 但仍然比Beta 控制圖的性能更好, 例如m= 20 時(shí), 誤差系數(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í),即參數(shù)β和α的先驗(yàn)值prior(β)和prior(α)小于其真實(shí)值時(shí),控制圖的性能隨著誤差系數(shù)的減小越來(lái)越好.

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

5 案例分析

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

表5 硬盤驅(qū)動(dòng)器失效時(shí)間Table 5 failure times of hard disk drives

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

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

6 結(jié)束語(yǔ)

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

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