譚 超，劉 堅，張 星

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

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

譚 超，劉 堅，張 星

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

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

貝葉斯均值控制圖；參數(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)濟性設計，并證明了其在經(jīng)濟性上的優(yōu)越性。朱慧明[10]研究了貝葉斯序貫過程質(zhì)量監(jiān)控模型，提出了基于多階段預報分布的貝葉斯多變量均值向量監(jiān)控模型。Marcellus[11，12]從統(tǒng)計控制角度，研究了均值漂移過程的貝葉斯分析問題。

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

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

1 貝葉斯均值控制圖

Marcellus[11，12]提出了貝葉斯均值控制圖以監(jiān)測均值的偏移。假設受控狀態(tài)下，獨立同分布的隨機樣本，其發(fā)生故障的時間間隔服從均值為的指數(shù)分布，則過程在T時間內(nèi)以概率=1-e-λHi從受控狀態(tài)進入到失控狀態(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)進入失控狀態(tài)，則可以得到在Hi-1到Hi過程中，從受控狀態(tài)進入第一類失控狀態(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）當π1(i)+π2(i)3q時，失控報警，檢查原因。并重置π0(0)=1、π1(0)=0和π2(0)=0。

（2）當π1(i)+π2(i)

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

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

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

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

當Phase I結(jié)束后，在phase II本文選取 μ0=0，σ0=1。λ=0.01。假設發(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，對應的μ2=-0.4，-0.6，-0.8，-1。

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

4 結(jié)論

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

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

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

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

（3）本文用“二分查找”設計了在參數(shù)估計條件下不同參數(shù)組合ARL0=200的控制線q，并用最小二乘估計，得出了一個關于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]朱慧明，管皓云，林靜等.基于多階段預報分布的貝葉斯多變量均值向量監(jiān)控模型[J].湖南大學學報（自然科學版），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).

（責任編輯/易永生）

O213.1

A

1002-6487（2016）21-0022-04