基于参数估计的贝叶斯均值控制图研究

谭 超，刘 坚，张 星

（湖南大学机械与运载工程学院，长沙 410082）

基于参数估计的贝叶斯均值控制图研究

谭 超，刘 坚，张 星

（湖南大学机械与运载工程学院，长沙 410082）

在传统的贝叶斯控制图研究中，通常将参数设为已知常量，忽略了工程实践性。文章对参数估计条件下的贝叶斯均值控制图性能展开研究，应用Monte Carlo仿真方法计算贝叶斯控制图的性能指标平均运行链长ARL0，分析了参数估计对控制图的性能影响，阐述了影响趋势的本质原因，设计了参数估计条件下保证统计性能的贝叶斯控制图设计方法，为拓展贝叶斯控制图的工程实践性提供了有益的尝试。

贝叶斯均值控制图；参数估计；平均运行链长

0 引言

贝叶斯控制图起源于Girshick和Rubin[1]、Bathers[2]和Taylor[3，4]的研究成果。Tagaras[5，6]提出了变参数的单边贝叶斯控制图，结合成本模型，提出了短周期运行过程的动态贝叶斯控制图。Makis[7]研究了短周期运行的最优多元贝叶斯控制问题。Nenes[8，9]研究了双边贝叶斯控制图的经济性设计，并证明了其在经济性上的优越性。朱慧明[10]研究了贝叶斯序贯过程质量监控模型，提出了基于多阶段预报分布的贝叶斯多变量均值向量监控模型。Marcellus[11，12]从统计控制角度，研究了均值漂移过程的贝叶斯分析问题。

上述有关贝叶斯控制图的研究均假定过程参数为已知。在工程实践过程中，过程参数通常是难以精确获取的，所以需要对过程参数进行参数估计。一般的控制图构建包含了两个阶段，Phase I和Phase II。在Phase I，主要目的是在受控状态下进行参数估计以构建控制图，而Phase II的主要目的是为了尽快检测出失控状态。

由于贝叶斯控制图的研究工作起步较晚，参数估计条件下的贝叶斯控制图性能研究工作鲜有报道，鉴于此，本文围绕参数估计条件下的贝叶斯控制图设计问题展开研究，以期回答如下三个问题：(1)估计参数替代已知参数对贝叶斯控制图的影响；(2)Phase I需要多大的样本值才能保证Phase II的性能；(3)Phase II的控制线应如何调整，以补偿Phase I样本值不够的情况。

1 贝叶斯均值控制图

Marcellus[11，12]提出了贝叶斯均值控制图以监测均值的偏移。假设受控状态下，独立同分布的随机样本，其发生故障的时间间隔服从均值为的指数分布，则过程在T时间内以概率=1-e-λHi从受控状态进入到失控状态。在T时间内，第一类失控状态为样本均值Yi向上偏移至μ1，其发生的概率为 p1，而第二类失控状态为样本均值Yi向下偏移至μ2，其发生的概率为 p2。过程在 Hi时刻抽样，i³1，Hi>Hi-1。在抽样时间Hi，抽取n个随机变量的均值Yi被检测得到。过程处于受控状态时，Yi的概率密度函数为。过程处于失控状态时，第一类失控状态和第二类失控状态的概率密度函数分别为

在初始时刻，π0(0)=1、π1(0)=0和π2(0)=0分别代表过程处于受控状态、第一类失控状态和第二类失控状态的概率。在Hi时刻，可以计算获得π0(i)、π1(i)和π2(i)，其分别是基于先验概率π0(i-1)、π1(i-1)和π2(i-1)的后验概率。同时可以获得两个新信息：过程又运行了Hi-Hi-1的时间，以及Yi。

假定过程从Hi-1时刻到Hi时刻由受控状态进入失控状态，则可以得到在Hi-1到Hi过程中，从受控状态进入第一类失控状态和第二类失控状态的概率分别是：

因此，在此过程中，处于受控状态的概率为：

在Hi时刻之前，处于第一类失控状态和第二类失控状态的概率分别为π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)³q时，失控报警，检查原因。并重置π0(0)=1、π1(0)=0和π2(0)=0。

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

2 参数估计对控制图性能的影响

当用估计参数取代已知参数时，贝叶斯均值控制图应考虑到参数估计量的变化性，即参数估计对控制图的影响。如果研究者没有考虑到可变化性，那么控制图的性能将会被严重的影响。估计量的精确性会随着Phase I所抽取的样本数的增加而增加。本文运用蒙特卡洛法仿真计算出贝叶斯均值控制图的受控时平均运行链长（In-Control Average Running Length，简称ARL0），以其表征参数估计对贝叶斯均值控制图性能的影响。

在以往的研究中，贝叶斯均值控制图是在μ0和σ0假设为已知的条件下设计的。当μ0和σ0未知时，由Albers等[13]的研究知，Phase I可以抽取m个独立的样本，用这m个样本的均值和样本标准差S分别对μ0和σ0进行估计。其表达式分别为：

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

当Phase I结束后，在phase II本文选取 μ0=0，σ0=1。λ=0.01。假设发生两类失控状态的可能性相同，那么p1=0.5、p2=0.5。由于不同的抽样时间间隔Hi和对贝叶斯均值控制图的影响相对较小[11，12]，因而本文的研究中Hi=1。影响贝叶斯均值控制图性能的主要参数是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）当参数设定为已知，且在n、μ1和μ2取不同值时，通过大样本的蒙特卡洛仿真计算出贝叶斯均值控制图的ARL0=200时对应的控制线q值，计算结果如表1所示。

（3）在Phase II，在第i个抽样时间间隔Hi，抽取n个样本，得到样本均值Yi。

（4）通过式（4）和式（5）分别更新π1(i)和π2(i)，并将π1(i)+π2(i)与步骤（1）得到的q进行比较。

（5）重复步骤（3）和步骤（4），直到发出失控的信号。记录一次运行链长。

（6）重复步骤（2）到步骤（5）50000次，得到ARL0。

表1 参数已知条件下ARL0=200的控制线q值

表2至表4给出了在参数估计的情况下，控制线为表1中所示q时，不同的m、n、μ1和μ2取值对应计算出来的ARL0值。根据表2至表4的结果表明，用参数已知情况下的控制线q作为参数估计情况下的控制线去仿真，贝叶斯均值控制图的性能会被严重的影响。如表3所示，当n=5，m=20，μ1=1和μ2=-1时，ARL0=374.1，这与参数已知情况下ARL0=200有87.5%的偏差。这个例子表明了，在建立控制图的过程中，不能忽视参数的可变性，否则会造成控制图性能的严重偏差。经对比表2至表4的结果，可以得出以下几个推论：

（1）随着 μ1与μ2之间的距离增大，ARL0会随之增长。如表3所示，当m=20时，随着μ1与μ2之间的距离增大，ARL0从205.3增大到374.1。这是因为在其他同等参数条件下，随着偏差允许范围增加，在μ0=0的条件下发生误报警的可能性降低，所以ARL0也会随之增加；

（2）随着n的增大，ARL0的大小会相应的增加。例如当 μ1=1，μ2=-1，m=20时，随着n的增大，ARL0从214.6增大到644.0。这是因为随着n的增加，样本均值会更加精确，因此误报的可能性会降低，ARL0会越来越大；

（3）随着m的增大，ARL0会随之减小，ARL0趋近参数已知情况。例如，表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时，m³40基本可以满足Phase I参数估计的样本需求；当 μ1=0.6，μ2=-0.6时m³100基本可以满足Phase I参数估计的样本需求；当μ1=0.8，μ2=-0.8时，m³400基本可以满足Phase I参数估计的样本需求；当μ1=1，μ2=-1时，m³1200基本可以满足Phase I参数估计的样本需求。由此可见，随着μ1与μ2之间的距离增大，Phase I所需要的m值也随之增加。因为推论（1），可得μ1与μ2之间的距离增大会导致ARL0的增大，从而需要更大的m，以保证估计的精确性。

表2 参数估计条件下n=1时的ARL0

表3 参数估计条件下n=5时的ARL0

表4 参数估计条件下n=10时的ARL0

3 参数估计条件的贝叶斯均值控制图控制线的设计

在某些实际应用过程中，可收集的实验样本相当多，因此在Phase I等待直到收集到1600或大于1600个样本是可行的。然而，在大多数的实际应用中，从经济性和实践性的角度来说，能够在Phase I收集到足够大的样本数是不可行的。因此，许多学者开始对质量控制图进行设计，使得质量控制图能够不需要假设参数。例如，Quesenberry[14]等的研究，这些研究的主要思想是对控制线进行设计，将控制线放宽，以达到满意的控制图的性能。然而，贝叶斯均值控制图的思路应与之相反，随着m减小时，应该将控制线q也相应的减小，因为m较小时，ARL0会相应的增加，所以，需要减小q，使得ARL0降低，以达到所需的ARL0。本文将根据“二分查找”，运用程序计算出在Phase I阶段，m、n、μ1和μ2为不同参数组合时，ARL0= 200时的控制线q，根据所得到的q，用最小二乘估计得出一个线性回归模型，以用来估计不同参数组合时所需要的最合适的q。

表5至表7给出了当n=1，5，10时，不同的m、μ1和μ2的参数组合满足ARL0=200需要的控制线q。这些控制线是由第二节所用的仿真步骤（2）至步骤（5）计算出来，唯一在其中增加了“二分查找”。“二分查找”的主要思想是在控制线q的一定范围内不断对q折半，直到寻找到一个q能够计算出满意的ARL0。

表5 参数估计条件下n=1时ARL0=200的控制线q值

表6 参数估计条件下n=5时ARL0=200的控制线q值

表7 参数估计条件下n=10时ARL0=200的控制线q值

表5至表7的结果表明改进的控制线的大小取决于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与参数已知时的q几乎相同，再次证明了Phase I的大样本能确保参数估计的精确性。然而，在实际情况中，m的值不同于表5至表7中给出的数值时，基于此，贝叶斯均值控制图的使用人员可以用插值法得到一个合适的改进控制线。而本文用表5至表7中给出的数据，对q进行最小二乘估计，计算得出一个形如a+blog10m的简单的线性回归模型，以计算不同m值所需要的改进控制线，如表8所示。例如，当n=5，μ1=0.6，μ2=-0.6时，改进控制线的线性回归函数为：

当m=150时，q=0.174，仿真可得ARL0=197.6，与ARL0=200之间仅有1.2%的误差。

4 结论

本文围绕参数估计条件下的贝叶斯均值控制图设计问题展开研究，通过样本的均值和样本标准差S分别对μ0和σ0进行参数估计，研究了参数估计对Marcellus[12，13]提出的贝叶斯均值控制图的影响。得出了如下结论：

表8 参数估计条件下ARL0=200的控制线q的线性回归模型

（1）在参数已知条件下，本文通过大样本的蒙特卡洛仿真计算出贝叶斯均值控制图不同参数组合下ARL0=200时对应的控制线q值，用其作为参数估计条件下对应参数组合的控制线。继而，在参数估计条件下进行仿真，结果发现ARL0与200之间产生了较大的偏差，证明了估计参数对贝叶斯均值控制图的性能确有较大的影响。

（2）本文在结论（1）的基础之上，给出了贝叶斯均值控制图在参数估计条件下不同参数组合在Phase I所需的m值。当μ1=0.4，μ2=-0.4时，m³40基本可以满足Phase I参数估计的样本需求；当μ1=0.6，μ2=-0.6时m³100基本可以满足Phase I参数估计的样本需求；当μ1=0.8，μ2=-0.8时，m³400基本可以满足Phase I参数估计的样本需求；当μ1=1，μ2=-1时，m³1200基本可以满足Phase I参数估计的样本需求。

（3）本文用“二分查找”设计了在参数估计条件下不同参数组合ARL0=200的控制线q，并用最小二乘估计，得出了一个关于q的线性回归函数，以补偿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]朱慧明，管皓云，林静等.基于多阶段预报分布的贝叶斯多变量均值向量监控模型[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