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

譚 超，劉 堅，張 星

（湖南大學(xué)機械與運載工程學(xué)院，長沙 410082）

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

譚 超，劉 堅，張 星

（湖南大學(xué)機械與運載工程學(xué)院，長沙 410082）

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

貝葉斯均值控制圖；參數(shù)估計；平均運行鏈長

0 引言

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

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

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

1 貝葉斯均值控制圖

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

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

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

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

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

其中：

在Hi抽樣后，計算出后驗π1(i)和π2(i)，選擇一個控制線q，0

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

（1）隨著 μ1與μ2之間的距離增大，ARL0會隨之增長。如表3所示，當(dāng)m=20時，隨著μ1與μ2之間的距離增大，ARL0從205.3增大到374.1。這是因為在其他同等參數(shù)條件下，隨著偏差允許范圍增加，在μ0=0的條件下發(fā)生誤報警的可能性降低，所以ARL0也會隨之增加；

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

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

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

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

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

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

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

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

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

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

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

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

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

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

4 結(jié)論

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

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

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

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

（3）本文用“二分查找”設(shè)計了在參數(shù)估計條件下不同參數(shù)組合ARL0=200的控制線q，并用最小二乘估計，得出了一個關(guān)于q的線性回歸函數(shù)，以補償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ù)報分布的貝葉斯多變量均值向量監(jiān)控模型[J].湖南大學(xué)學(xué)報（自然科學(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