譚 超，劉 堅(jiān)，張 星

（湖南大學(xué)機(jī)械與運(yùn)載工程學(xué)院，長(zhǎng)沙 410082）

基于參數(shù)估計(jì)的貝葉斯均值控制圖研究

譚 超，劉 堅(jiān)，張 星

（湖南大學(xué)機(jī)械與運(yùn)載工程學(xué)院，長(zhǎng)沙 410082）

在傳統(tǒng)的貝葉斯控制圖研究中，通常將參數(shù)設(shè)為已知常量，忽略了工程實(shí)踐性。文章對(duì)參數(shù)估計(jì)條件下的貝葉斯均值控制圖性能展開(kāi)研究，應(yīng)用Monte Carlo仿真方法計(jì)算貝葉斯控制圖的性能指標(biāo)平均運(yùn)行鏈長(zhǎng)ARL0，分析了參數(shù)估計(jì)對(duì)控制圖的性能影響，闡述了影響趨勢(shì)的本質(zhì)原因，設(shè)計(jì)了參數(shù)估計(jì)條件下保證統(tǒng)計(jì)性能的貝葉斯控制圖設(shè)計(jì)方法，為拓展貝葉斯控制圖的工程實(shí)踐性提供了有益的嘗試。

貝葉斯均值控制圖；參數(shù)估計(jì)；平均運(yùn)行鏈長(zhǎng)

0 引言

貝葉斯控制圖起源于Girshick和Rubin[1]、Bathers[2]和Taylor[3，4]的研究成果。Tagaras[5，6]提出了變參數(shù)的單邊貝葉斯控制圖，結(jié)合成本模型，提出了短周期運(yùn)行過(guò)程的動(dòng)態(tài)貝葉斯控制圖。Makis[7]研究了短周期運(yùn)行的最優(yōu)多元貝葉斯控制問(wèn)題。Nenes[8，9]研究了雙邊貝葉斯控制圖的經(jīng)濟(jì)性設(shè)計(jì)，并證明了其在經(jīng)濟(jì)性上的優(yōu)越性。朱慧明[10]研究了貝葉斯序貫過(guò)程質(zhì)量監(jiān)控模型，提出了基于多階段預(yù)報(bào)分布的貝葉斯多變量均值向量監(jiān)控模型。Marcellus[11，12]從統(tǒng)計(jì)控制角度，研究了均值漂移過(guò)程的貝葉斯分析問(wèn)題。

上述有關(guān)貝葉斯控制圖的研究均假定過(guò)程參數(shù)為已知。在工程實(shí)踐過(guò)程中，過(guò)程參數(shù)通常是難以精確獲取的，所以需要對(duì)過(guò)程參數(shù)進(jìn)行參數(shù)估計(jì)。一般的控制圖構(gòu)建包含了兩個(gè)階段，Phase I和Phase II。在Phase I，主要目的是在受控狀態(tài)下進(jìn)行參數(shù)估計(jì)以構(gòu)建控制圖，而Phase II的主要目的是為了盡快檢測(cè)出失控狀態(tài)。

由于貝葉斯控制圖的研究工作起步較晚，參數(shù)估計(jì)條件下的貝葉斯控制圖性能研究工作鮮有報(bào)道，鑒于此，本文圍繞參數(shù)估計(jì)條件下的貝葉斯控制圖設(shè)計(jì)問(wèn)題展開(kāi)研究，以期回答如下三個(gè)問(wèn)題：(1)估計(jì)參數(shù)替代已知參數(shù)對(duì)貝葉斯控制圖的影響；(2)Phase I需要多大的樣本值才能保證Phase II的性能；(3)Phase II的控制線(xiàn)應(yīng)如何調(diào)整，以補(bǔ)償Phase I樣本值不夠的情況。

1 貝葉斯均值控制圖

Marcellus[11，12]提出了貝葉斯均值控制圖以監(jiān)測(cè)均值的偏移。假設(shè)受控狀態(tài)下，獨(dú)立同分布的隨機(jī)樣本，其發(fā)生故障的時(shí)間間隔服從均值為的指數(shù)分布，則過(guò)程在T時(shí)間內(nèi)以概率=1-e-λHi從受控狀態(tài)進(jìn)入到失控狀態(tài)。在T時(shí)間內(nèi)，第一類(lèi)失控狀態(tài)為樣本均值Yi向上偏移至μ1，其發(fā)生的概率為 p1，而第二類(lèi)失控狀態(tài)為樣本均值Yi向下偏移至μ2，其發(fā)生的概率為 p2。過(guò)程在 Hi時(shí)刻抽樣，i31，Hi>Hi-1。在抽樣時(shí)間Hi，抽取n個(gè)隨機(jī)變量的均值Yi被檢測(cè)得到。過(guò)程處于受控狀態(tài)時(shí)，Yi的概率密度函數(shù)為。過(guò)程處于失控狀態(tài)時(shí)，第一類(lèi)失控狀態(tài)和第二類(lèi)失控狀態(tài)的概率密度函數(shù)分別為

在初始時(shí)刻，π0(0)=1、π1(0)=0和π2(0)=0分別代表過(guò)程處于受控狀態(tài)、第一類(lèi)失控狀態(tài)和第二類(lèi)失控狀態(tài)的概率。在Hi時(shí)刻，可以計(jì)算獲得π0(i)、π1(i)和π2(i)，其分別是基于先驗(yàn)概率π0(i-1)、π1(i-1)和π2(i-1)的后驗(yàn)概率。同時(shí)可以獲得兩個(gè)新信息：過(guò)程又運(yùn)行了Hi-Hi-1的時(shí)間，以及Yi。

假定過(guò)程從Hi-1時(shí)刻到Hi時(shí)刻由受控狀態(tài)進(jìn)入失控狀態(tài)，則可以得到在Hi-1到Hi過(guò)程中，從受控狀態(tài)進(jìn)入第一類(lèi)失控狀態(tài)和第二類(lèi)失控狀態(tài)的概率分別是：

因此，在此過(guò)程中，處于受控狀態(tài)的概率為：

在Hi時(shí)刻之前，處于第一類(lèi)失控狀態(tài)和第二類(lèi)失控狀態(tài)的概率分別為π0(i-1)a1(i)+π1(i-1)和π0(i-1)a2(i)+ π2(i-1)。因此，在抽取了Yi后，后驗(yàn)失控概率為：

其中：

在Hi抽樣后，計(jì)算出后驗(yàn)π1(i)和π2(i)，選擇一個(gè)控制線(xiàn)q，0

（1）當(dāng)π1(i)+π2(i)3q時(shí)，失控報(bào)警，檢查原因。并重置π0(0)=1、π1(0)=0和π2(0)=0。

（2）當(dāng)π1(i)+π2(i)

2 參數(shù)估計(jì)對(duì)控制圖性能的影響

當(dāng)用估計(jì)參數(shù)取代已知參數(shù)時(shí)，貝葉斯均值控制圖應(yīng)考慮到參數(shù)估計(jì)量的變化性，即參數(shù)估計(jì)對(duì)控制圖的影響。如果研究者沒(méi)有考慮到可變化性，那么控制圖的性能將會(huì)被嚴(yán)重的影響。估計(jì)量的精確性會(huì)隨著Phase I所抽取的樣本數(shù)的增加而增加。本文運(yùn)用蒙特卡洛法仿真計(jì)算出貝葉斯均值控制圖的受控時(shí)平均運(yùn)行鏈長(zhǎng)（In-Control Average Running Length，簡(jiǎn)稱(chēng)ARL0），以其表征參數(shù)估計(jì)對(duì)貝葉斯均值控制圖性能的影響。

在以往的研究中，貝葉斯均值控制圖是在μ0和σ0假設(shè)為已知的條件下設(shè)計(jì)的。當(dāng)μ0和σ0未知時(shí)，由Albers等[13]的研究知，Phase I可以抽取m個(gè)獨(dú)立的樣本，用這m個(gè)樣本的均值和樣本標(biāo)準(zhǔn)差S分別對(duì)μ0和σ0進(jìn)行估計(jì)。其表達(dá)式分別為：

由 Xi～N(μ0，σ02)，可 知由，那么可推得，S～即

當(dāng)Phase I結(jié)束后，在phase II本文選取 μ0=0，σ0=1。λ=0.01。假設(shè)發(fā)生兩類(lèi)失控狀態(tài)的可能性相同，那么p1=0.5、p2=0.5。由于不同的抽樣時(shí)間間隔Hi和對(duì)貝葉斯均值控制圖的影響相對(duì)較小[11，12]，因而本文的研究中Hi=1。影響貝葉斯均值控制圖性能的主要參數(shù)是n、μ1和μ2。本文選取n=1，5，10，μ1=0.4，0.6，0.8，1，對(duì)應(yīng)的μ2=-0.4，-0.6，-0.8，-1。

本文研究不同的m、n、μ1和μ2對(duì)貝葉斯均值控制圖ARL0的影響，其研究步驟如下：

（1）當(dāng)參數(shù)設(shè)定為已知，且在n、μ1和μ2取不同值時(shí)，通過(guò)大樣本的蒙特卡洛仿真計(jì)算出貝葉斯均值控制圖的ARL0=200時(shí)對(duì)應(yīng)的控制線(xiàn)q值，計(jì)算結(jié)果如表1所示。

（3）在Phase II，在第i個(gè)抽樣時(shí)間間隔Hi，抽取n個(gè)樣本，得到樣本均值Yi。

（4）通過(guò)式（4）和式（5）分別更新π1(i)和π2(i)，并將π1(i)+π2(i)與步驟（1）得到的q進(jìn)行比較。

（5）重復(fù)步驟（3）和步驟（4），直到發(fā)出失控的信號(hào)。記錄一次運(yùn)行鏈長(zhǎng)。

（6）重復(fù)步驟（2）到步驟（5）50000次，得到ARL0。

表1 參數(shù)已知條件下ARL0=200的控制線(xiàn)q值

表2至表4給出了在參數(shù)估計(jì)的情況下，控制線(xiàn)為表1中所示q時(shí)，不同的m、n、μ1和μ2取值對(duì)應(yīng)計(jì)算出來(lái)的ARL0值。根據(jù)表2至表4的結(jié)果表明，用參數(shù)已知情況下的控制線(xiàn)q作為參數(shù)估計(jì)情況下的控制線(xiàn)去仿真，貝葉斯均值控制圖的性能會(huì)被嚴(yán)重的影響。如表3所示，當(dāng)n=5，m=20，μ1=1和μ2=-1時(shí)，ARL0=374.1，這與參數(shù)已知情況下ARL0=200有87.5%的偏差。這個(gè)例子表明了，在建立控制圖的過(guò)程中，不能忽視參數(shù)的可變性，否則會(huì)造成控制圖性能的嚴(yán)重偏差。經(jīng)對(duì)比表2至表4的結(jié)果，可以得出以下幾個(gè)推論：

（1）隨著 μ1與μ2之間的距離增大，ARL0會(huì)隨之增長(zhǎng)。如表3所示，當(dāng)m=20時(shí)，隨著μ1與μ2之間的距離增大，ARL0從205.3增大到374.1。這是因?yàn)樵谄渌葏?shù)條件下，隨著偏差允許范圍增加，在μ0=0的條件下發(fā)生誤報(bào)警的可能性降低，所以ARL0也會(huì)隨之增加；

（2）隨著n的增大，ARL0的大小會(huì)相應(yīng)的增加。例如當(dāng) μ1=1，μ2=-1，m=20時(shí)，隨著n的增大，ARL0從214.6增大到644.0。這是因?yàn)殡S著n的增加，樣本均值會(huì)更加精確，因此誤報(bào)的可能性會(huì)降低，ARL0會(huì)越來(lái)越大；

（3）隨著m的增大，ARL0會(huì)隨之減小，ARL0趨近參數(shù)已知情況。例如，表2所示，當(dāng)μ1=0.6，μ2=-0.6時(shí)，隨著m的增加，ARL0逐漸趨近200。因?yàn)楫?dāng)m取值小時(shí)，μ0和σ0會(huì)因估計(jì)不精確導(dǎo)致波動(dòng)較大，從而ARL0的取值會(huì)呈1：1的比例在200上下波動(dòng)，但因大于200的波動(dòng)幅度遠(yuǎn)大于小于200的幅度，從而取均值時(shí)，ARL0會(huì)大于200。而當(dāng)m取值越來(lái)越大，μ0和σ0的估計(jì)會(huì)越來(lái)越精確，使μ0和σ0的波動(dòng)越來(lái)越小，從而使得ARL0逐漸趨近200；

（4）當(dāng) μ1=0.4，μ2=-0.4時(shí)，m340基本可以滿(mǎn)足Phase I參數(shù)估計(jì)的樣本需求；當(dāng) μ1=0.6，μ2=-0.6時(shí)m3100基本可以滿(mǎn)足Phase I參數(shù)估計(jì)的樣本需求；當(dāng)μ1=0.8，μ2=-0.8時(shí)，m3400基本可以滿(mǎn)足Phase I參數(shù)估計(jì)的樣本需求；當(dāng)μ1=1，μ2=-1時(shí)，m31200基本可以滿(mǎn)足Phase I參數(shù)估計(jì)的樣本需求。由此可見(jiàn)，隨著μ1與μ2之間的距離增大，Phase I所需要的m值也隨之增加。因?yàn)橥普摚?），可得μ1與μ2之間的距離增大會(huì)導(dǎo)致ARL0的增大，從而需要更大的m，以保證估計(jì)的精確性。

表2 參數(shù)估計(jì)條件下n=1時(shí)的ARL0

表3 參數(shù)估計(jì)條件下n=5時(shí)的ARL0

表4 參數(shù)估計(jì)條件下n=10時(shí)的ARL0

3 參數(shù)估計(jì)條件的貝葉斯均值控制圖控制線(xiàn)的設(shè)計(jì)

在某些實(shí)際應(yīng)用過(guò)程中，可收集的實(shí)驗(yàn)樣本相當(dāng)多，因此在Phase I等待直到收集到1600或大于1600個(gè)樣本是可行的。然而，在大多數(shù)的實(shí)際應(yīng)用中，從經(jīng)濟(jì)性和實(shí)踐性的角度來(lái)說(shuō)，能夠在Phase I收集到足夠大的樣本數(shù)是不可行的。因此，許多學(xué)者開(kāi)始對(duì)質(zhì)量控制圖進(jìn)行設(shè)計(jì)，使得質(zhì)量控制圖能夠不需要假設(shè)參數(shù)。例如，Quesenberry[14]等的研究，這些研究的主要思想是對(duì)控制線(xiàn)進(jìn)行設(shè)計(jì)，將控制線(xiàn)放寬，以達(dá)到滿(mǎn)意的控制圖的性能。然而，貝葉斯均值控制圖的思路應(yīng)與之相反，隨著m減小時(shí)，應(yīng)該將控制線(xiàn)q也相應(yīng)的減小，因?yàn)閙較小時(shí)，ARL0會(huì)相應(yīng)的增加，所以，需要減小q，使得ARL0降低，以達(dá)到所需的ARL0。本文將根據(jù)“二分查找”，運(yùn)用程序計(jì)算出在Phase I階段，m、n、μ1和μ2為不同參數(shù)組合時(shí)，ARL0= 200時(shí)的控制線(xiàn)q，根據(jù)所得到的q，用最小二乘估計(jì)得出一個(gè)線(xiàn)性回歸模型，以用來(lái)估計(jì)不同參數(shù)組合時(shí)所需要的最合適的q。

表5至表7給出了當(dāng)n=1，5，10時(shí)，不同的m、μ1和μ2的參數(shù)組合滿(mǎn)足ARL0=200需要的控制線(xiàn)q。這些控制線(xiàn)是由第二節(jié)所用的仿真步驟（2）至步驟（5）計(jì)算出來(lái)，唯一在其中增加了“二分查找”?！岸植檎摇钡闹饕枷胧窃诳刂凭€(xiàn)q的一定范圍內(nèi)不斷對(duì)q折半，直到尋找到一個(gè)q能夠計(jì)算出滿(mǎn)意的ARL0。

表5 參數(shù)估計(jì)條件下n=1時(shí)ARL0=200的控制線(xiàn)q值

表6 參數(shù)估計(jì)條件下n=5時(shí)ARL0=200的控制線(xiàn)q值

表7 參數(shù)估計(jì)條件下n=10時(shí)ARL0=200的控制線(xiàn)q值

表5至表7的結(jié)果表明改進(jìn)的控制線(xiàn)的大小取決于m、n、μ1和μ2。較小的m需要較小的控制線(xiàn)q去得到滿(mǎn)意的ARL0。例如表5所示，當(dāng)n=10，μ1=0.6，μ2=-0.6時(shí)，從m=20至m=1600的過(guò)程中ARL0從0.100增加到0.132。且當(dāng)m=1600時(shí)，q與參數(shù)已知時(shí)的q幾乎相同，再次證明了Phase I的大樣本能確保參數(shù)估計(jì)的精確性。然而，在實(shí)際情況中，m的值不同于表5至表7中給出的數(shù)值時(shí)，基于此，貝葉斯均值控制圖的使用人員可以用插值法得到一個(gè)合適的改進(jìn)控制線(xiàn)。而本文用表5至表7中給出的數(shù)據(jù)，對(duì)q進(jìn)行最小二乘估計(jì)，計(jì)算得出一個(gè)形如a+blog10m的簡(jiǎn)單的線(xiàn)性回歸模型，以計(jì)算不同m值所需要的改進(jìn)控制線(xiàn)，如表8所示。例如，當(dāng)n=5，μ1=0.6，μ2=-0.6時(shí)，改進(jìn)控制線(xiàn)的線(xiàn)性回歸函數(shù)為：

當(dāng)m=150時(shí)，q=0.174，仿真可得ARL0=197.6，與ARL0=200之間僅有1.2%的誤差。

4 結(jié)論

本文圍繞參數(shù)估計(jì)條件下的貝葉斯均值控制圖設(shè)計(jì)問(wèn)題展開(kāi)研究，通過(guò)樣本的均值和樣本標(biāo)準(zhǔn)差S分別對(duì)μ0和σ0進(jìn)行參數(shù)估計(jì)，研究了參數(shù)估計(jì)對(duì)Marcellus[12，13]提出的貝葉斯均值控制圖的影響。得出了如下結(jié)論：

表8 參數(shù)估計(jì)條件下ARL0=200的控制線(xiàn)q的線(xiàn)性回歸模型

（1）在參數(shù)已知條件下，本文通過(guò)大樣本的蒙特卡洛仿真計(jì)算出貝葉斯均值控制圖不同參數(shù)組合下ARL0=200時(shí)對(duì)應(yīng)的控制線(xiàn)q值，用其作為參數(shù)估計(jì)條件下對(duì)應(yīng)參數(shù)組合的控制線(xiàn)。繼而，在參數(shù)估計(jì)條件下進(jìn)行仿真，結(jié)果發(fā)現(xiàn)ARL0與200之間產(chǎn)生了較大的偏差，證明了估計(jì)參數(shù)對(duì)貝葉斯均值控制圖的性能確有較大的影響。

（2）本文在結(jié)論（1）的基礎(chǔ)之上，給出了貝葉斯均值控制圖在參數(shù)估計(jì)條件下不同參數(shù)組合在Phase I所需的m值。當(dāng)μ1=0.4，μ2=-0.4時(shí)，m340基本可以滿(mǎn)足Phase I參數(shù)估計(jì)的樣本需求；當(dāng)μ1=0.6，μ2=-0.6時(shí)m3100基本可以滿(mǎn)足Phase I參數(shù)估計(jì)的樣本需求；當(dāng)μ1=0.8，μ2=-0.8時(shí)，m3400基本可以滿(mǎn)足Phase I參數(shù)估計(jì)的樣本需求；當(dāng)μ1=1，μ2=-1時(shí)，m31200基本可以滿(mǎn)足Phase I參數(shù)估計(jì)的樣本需求。

（3）本文用“二分查找”設(shè)計(jì)了在參數(shù)估計(jì)條件下不同參數(shù)組合ARL0=200的控制線(xiàn)q，并用最小二乘估計(jì)，得出了一個(gè)關(guān)于q的線(xiàn)性回歸函數(shù)，以補(bǔ)償Phase I樣本值m不夠的情況。

[1]Girshick M A，Rubin H.A Bayes Approach to a Quality Control Model [J].Annals of Mathematical Statistics，1952，23(1).

[2]Bather J A.Control Charts and Minimization of Costs[J].Journal of the Royal Statistical Society-Series B，1963，25(1).

[3]Taylor H M.Markovian Sequential Replacement Processes[J].Annals of Mathematical Statistics，1965，36(6).

[4]Taylor H M.Statistical Control of a Gaussian Process[J].Technomet? rics，1967，9(1).

[5]Tagaras G.A Dynamic Programming Approach to the Economic De?sign of-Charts[J].IIE Transactions，1994，26(3).

[6]Tagaras G.Dynamic Control Charts for Finite Production Runs[J].Eu?ropean Journal of Operational Research，1996，91(1).

[7]Makis V.Multivariate Bayesian Process Control for a Finite Produc?tion Run[J].European Journal of Operational Research，2009，194(3).

[8]Nenes G，Tagaras G.The Economically Designed Two-Sided Bayes?ianControl Charts[J].European Journal of Operational Research， 2007，183(1).

[9]Nenes G.A New Approach for the Economic Design of Fully Adaptive Control Charts[J].International Journal of Production Economics， 2011，131(2).

[10]朱慧明，管皓云，林靜等.基于多階段預(yù)報(bào)分布的貝葉斯多變量均值向量監(jiān)控模型[J].湖南大學(xué)學(xué)報(bào)（自然科學(xué)版），2011，38(3).

[11]Marcellus R L.Bayesian Statistical Process Control[J].Quality Engi?neering，2007，20(1).

[12]Marcellus R L.Bayesian Monitoring to Detect a Shift in Process Mean[J].Quality and Reliability Engineering International，2008，24 (3).

[13]Albers W，Kallenberg W C M.Are Estimated Control Charts in Con?trol[J].Statistics，2004，38(1).

[14]Quesenberry C P.The Effect of Sample Size on Estimated Limits forand X Control Chart[J].Journal of Quality Technology，1993，25 (4).

（責(zé)任編輯/易永生）

O213.1

A

1002-6487（2016）21-0022-04