基于最大熵分布的控制图改进与评价研究

宋明顺，杨 铭，方兴华

(中国计量大学经济与管理学院，浙江 杭州 310018)

1 引言

以休哈特(Shewhart)控制图为代表的控制图技术在统计过程控制(Statistical Process Control，SPC)中具有很好的应用，许多研究已经证明Shewhart控制图对服从正态分布的质量特性参数具有较好的监控能力[1-2]。但在很多情况下，过程的质量特性参数并不服从正态分布，若此时仍按照正态分布构建控制图就会引起监控能力的下降。因此，研究一般分布情况下的控制图具有重要意义。

在控制图的研究中，多数学者主要是从控制图的构建和控制图效果的评价两方面展开。在控制图的构建方面，Page[3]基于正态分布，提出了累积和(Cumulative Sum，CUSUM)控制图，通过把每个样本值的偏移进行累加，可以更快的发现偏移。Roberts[4]在Page[3]的基础上，提出了指数加权移动平均(Exponentially Weighted Moving Average，EWMA)控制图，利用历史样本数据，同样可以监控过程发生的小偏移。Tran等[5]基于CUSUM控制图对二元正态分布下的人口比例进行监控，证明了在小偏移的情况下CUSUM控制图具有更好的监控效果。Park等[6]提出了最小水平集合的理论，在给定第Ⅰ类错误α的情况下，依据最小水平集确定控制界限。在此基础上，Das和Zhou[7]、Fang Xinghua等[8]基于最大熵，将该理论应用于控制图的构建，与传统控制图相比，新构建的控制图可以不用假定分布的形式，更适合于一般分布的情况。

在控制图效果的评价方面，平均运行长度(Average Run Length, ARL)理解容易、计算简单，能够较好评价控制图的效果，因此作为控制图常用的评价指标。王兆军等[9]详细阐述了不同控制图ARL的计算方法，Sanusi等[10]联合Shewhart和CUSUM控制图构建了一种联合控制图，通过ARL对该控制图进行评价，给出了联合控制图在不同偏移情况下的应用。Shu Lianjie等[11]用马尔科夫链方法推导了自适应CUSUM控制图的ARL的计算公式，推广了马尔科夫链的应用范围。Vargas等[12]对比了CUSUM控制图和EWMA控制图，以ARL作为检测失控状态的指标，得出了两种控制图在不同控制区域下的应用情况。

综上所述，尽管有很多的学者在控制图的构建和评价方面进行了研究，但这些研究大部分建立在质量特性参数服从正态分布的基础上的。而当这种假设不满足时，应用这些控制图进行过程监测容易产生较大的偏差。因此，本研究提出一种建立在质量特性参数“分布自由型”基础上的控制图的构建方法。由于最大熵作为一种分布估计方法，不对分布进行人为假定，可以更好地反映质量特性参数服从的真实分布[13]，在借鉴Barzdajn[14]、Alwan等[15]、Chiang等[16]的研究之后，本文应用最大熵分布对质量特性参数进行拟合，并依此构建控制图和评价方法。

本文应用最大熵分布对控制图进行改进的基本思想为：首先确定质量特性参数的最大熵分布，然后根据最大熵分布对Shewhart控制图进行改进，并结合ARL作为评价指标，对改进后的控制图进行评价。最后，还提出了基于最大熵的CUSUM控制图的性能评价方法。本研究通过一个仿真案例来阐释最大熵对Shewhart控制图和CUSUM控制图的改进过程，以及对改进结果的比较。

2 基于最大熵分布的控制图的构建和评价

本节首先介绍了最大熵，然后提出基于最大熵的Shewhart控制图的构建方法和CUSUM控制图的评价方法。

2.1 最大熵原理

信息熵(Entropy)最早由Shannon[17]提出，用来表示信息的不确定性程度，熵越大，信息的不确定性程度越高，之后，Jaynes[18]在这基础上提出了最大熵(Maximum Entropy, ME)的概念，并且证明“在部分信息的基础上，熵值达到最大时的分布是随机变量唯一的无偏分布估计”，学者把基于以上原则确定分布的思想称为最大熵原理(Principle of Maximum Entropy, PME)。

近年来，应用PME确定随机变量分布已经广泛的应用于工程设计[19]、可靠性评估[20]、金融产品定价[21]等。通常，根据PME估计随机变量X的最大熵分布的过程可以转化为求解如下数学规划的问题：

(1)

其中，H(x)为变量X的熵函数，f(x)为变量X的最大熵分布概率密度函数，ui(x)表示X的若干函数，Mi为常数。约束条件一般从X已有的先验信息提取，并且进一步转化为数学函数ui(x)，而ui(x)的形式在很大程度上决定了公式(1)中规划模型的复杂程度，如ui(x)为原点矩、中心矩函数时，高阶矩时模型求解的复杂性要明显高于低阶矩的情况。

为了获得f(x)的表达式，引入拉格朗日乘子法[22]，得到f(x)为：

(2)

其中λ0,λ1,λ2,…,λm为拉格朗日乘子，为待定系数。为了进一步得到f(x)的显性表达式，本研究利用Matlab中的非线性规划功能求解λ0,λ1,λ2,…,λm。

2.2 最大熵分布对Shewhart控制图的改进

在控制图的构建中，控制界限是最重要的参数之一，图1给出了质量特性参数的控制界限示意图，其中UCL为上控制限，LCL为下控制限。传统的Shewhart控制图采用[μ-3σ,μ+3σ]作为最佳控制区域，其中μ和σ分别为质量特性参数的均值和标准差，但这个控制区域适用于质量特性值为正态分布的情况，而当质量特性参数偏移正态分布时，应用该控制区域进行过程监控时就会产生较大的风险。

图1 质量特性参数X控制区域构建的示意图

2.2.1 控制图的构建原则

Shewhart控制图构建的核心思想是当质量特性参数服从正态分布时，控制区域[μ-3σ,μ+3σ]的两类错误之和(第Ⅰ类错误和第Ⅱ类错误)最小，在实际的产品质量控制中，应用这个控制界限可以使“虚发报警”和“漏发报警”的概率之和达到最小，因而企业在质量控制中承担的风险最小，是一种“经济性”的做法。

当质量特性参数的分布为单峰分布时，并且应用最大熵分布对其进行拟合，遵循控制图的“经济性”原则，并借鉴Park等[6]的研究，将质量特性参数的监测分布总体分为两个部分：受控区域和失控区域。如果监测的样本在受控区域内，则称为过程受控，否则为过程失控(如图1所示)。其中受控区域A可表示为：

A={X|P(X∈A)≥1-α,X～F0}

(3)

其中α为第Ⅰ类错误概率(在图1中，在对称分布中α1=α2=0.5α)，X为用于过程监测的质量特性参数，通常为产品的技术特性参数，F0表示X的分布，概率密度函数为f(x)。在特定的α值下，从统计上可以找到无数多个满足公式(3)的集合A。根据Shewhart控制图的“经济性”原则，此时必须找到使第Ⅱ类错误β达到最小时的集合，我们称之为最佳受控区域A*，此时两类错误之和α+β达到最小。

当过程失控时，此时从统计学上看，X偏离了原分布F0，若监测的样本落在受控区域内，则会产生第Ⅱ类错误β，即“漏发报警”。假设失控状态下的分布为F1，相应的概率密度函数为f'(x)，则第Ⅱ类错误为：

(4)

在实际的产品质量监测的过程中，过程失控的原因往往是比较复杂的，可能由多种原因引起，因而f'(x)的表达式通常是未知的。为了便于计算，我们采用核密度估计对未知分布进行拟合。由于核密度估计被证明与原分布的一致性，而且能够较好的体现数据特征，具有较好的拟合效果[23]，这也是目前在分布估计中常用的一种方法。

(5)

其中K(·)表示核函数(非负，积分为1，均值为0，具有概率密度的性质)，hn为带宽，n为样本总量，Xi为单个质量特性参数的观测值。因此将公式(5)代入公式(4)，得到：

(6)

由于核密度估计时积分计算很复杂，为了简化求解，我们经常采用经验分布估计核密度积分计算[25]，其中经验分布函数Ln(A)为：

(7)

(8)

因此当质量特性参数为服从单峰分布的随机变量时，给定第Ⅰ类错误α，则当控制控制区间长度(UCL-LCL)取最短时，则两类错误之和α+β最小。

2.2.2 基于最大熵分布控制图的构建步骤

应用最大熵分布对Shewhart控制图进行改进，其改进的思想是用最大熵分布近似质量特性参数的分布，并根据“经济性”原则构建控制图，其主要的步骤如下：

(1)应用最大熵分布拟合质量特性参数X分布的概率密度函数f(X)，如公式(1)所示；

(2)给定第Ⅰ类错误α，寻找最佳控制区域A*=[LCL,UCL]以满足：

(9)

可以应用数值法等近似求解，并利用Matlab等软件实现；

(3)对过程进行监控，当观测样本落在受控区域外，则过程失控报警。

2.2.3 对改进Shewhart控制图的评价

在控制图效果的评价中，通常采用ARL作为评价标准。当过程受控时，ARL表示开始检测到发生报警信号平均抽取的样本数，此时我们希望其越大越好；反之，当过程失控时，ARL表示开始检测到发现失控信号平均抽取的样本数，此时其值越小越好。

对于Shewhart控制图，只需要知道第Ⅰ类错误α和第Ⅱ类错误β，就可以计算出受控状态下和失控状态下的ARL，如下所示：

(10)

(11)

其中ARL0表示过程受控时的平均运行长度，ARL1表示过程失控时的平均运行长度。由于本研究中第Ⅰ类错误α是一个确定值，对不同控制图没有影响，因此本研究应用ARL1来评价控制图的性能。

2.3 最大熵分布对CUSUM控制图评价方法的改进

一般来说，Shewhart控制图对过程大偏移的监控比较有效，而对于过程中的小偏移的监测不敏感，为此Page[3]、Tran等[5]基于序贯比检验，构建了更适用于小偏移监测的累积和(CUSUM)控制图。针对CUSUM控制图的评价，本部分提出了基于最大熵分布的控制图评价方法，并以自适应CUSUM控制图为例进行说明。

2.3.1 自适应CUSUM控制图

传统CUSUM控制图的原理是对每一个观测值的偏移进行累加，从而达到局部放大的效果。但在现实的生产中，很难获取偏移的确切值，这给控制图的设计带来了困扰，Shu Lianjie和Jiang Wei[26]提出了自适应CUSUM控制图，可以很好地解决这个问题。对于检测向上漂移，采用检测统计量为：

(12)

在上式中，ARL0表示受控状态下的平均运行长度，在试验中可取100-1000，本文中取ARL0=370。大量文献表明(如Sparks和Ross[27]，Huang等[28])，当k=δ/2，可以得到一个最优的CUSUM控制图。

2.3.2 基于最大熵分布的CUSUM控制图的评价

针对CUSUM控制图的评价，在计算ARL时，通常采用积分法、蒙特卡洛仿真法和马尔科夫链法[29]。相较于其他两种方法，马尔科夫链法适用于复杂的计算，同时准确度较高，因此本文采用马尔科夫链法计算CUSUM控制图的ARL。该方法的核心思想是把监控的质量特性参数近似看作一个状态有限的马尔可夫链，把质量特性参数的各取值区间作为马尔可夫链的各状态空间，以此构建转移概率矩阵，从而根据转移概率矩阵确定平均运行长度ARL的大小[30]。

对于一个马尔科夫链模型，主要由S、P两种元素构成，其中：

(1)S表示整个系统所有可能的状态所组成的非空的状态集，也可称为系统的状态空间，它可以是有限集合、可列集合或是有限可列集合。本文假定S是数集(即有限或可列)，用小写字母a,b或(Sa,Sb)来表示状态。

(2)P为系统的转移概率矩阵，Pab表示在时刻l为状态a，而在时刻l+1时为状态b的概率。

本文以向上的单侧偏移为例，假设受控区域为(0,H)，将该区域分为m个子区间，则每个子区间的宽度d=H/m，其中状态0到m-1为受控状态，状态m为失控状态，如图2所示。

图2 马尔科夫链状态转移空间

该马尔科夫链的一步转移概率矩阵P可表示为：

(13)

其中I为单位矩阵，1为元素全为1的列向量，R为m×m阶矩阵，R中元素表示由状态a转移到状态b的概率pab，即pab代表Ct-1落在状态a时，Ct落在状态b的概率，所以有：

(1)当b=0时，

pab=Pr(Ct∈S0|Ct-1∈Sa)=

Pr(Xa≤m-ad+d/2)=

F(Xa≤m-ad+d/2)

其中：a=1,2,3…,m-1，b=0

(2)当1

pab=Pr(Ct∈Sb|Ct-1∈Sa)=Pr[(b-a)d-d/2

(3)当b=m时，

pab=Pr(Ct∈Sb|Ct-1∈Sa)=Pr(X-K≤(m-a)d+d/2)=Pr(X≤(m-a+0.5)d+K)=F(X≤(m-a+0.5)d+K)

其中：a,b=1,2,3,…,m-1；b=m，F(·)为分布函数。

根据马尔科夫链的性质[31]，可以得到平均运行长度ARL的公式：

ARL=P0(I-R)-11

(14)

其中，P0是初始状态概率矩阵，为简单起见，设初始概率由中间点出发，即P0=[0,0,…,1,…,0,0]1×m。根据m的奇偶性，中间值的位置可分为：

(1)若m为偶数，则令P0为第m/2列为中间值，即P0(1,m/2)=1，其余值为0的向量。

(2)若m为奇数，则P0为中间值为1，其余值为0的向量。

由以上分析可以看出，在用马尔科夫链法计算CUSUM控制图的ARL时，需要知道观测值Xt的分布，且在大部分的研究中，ARL的计算建立在受控时过程参数都已知的假设基础上[32]，如假设均值μ0=0和标准差σ0=1的正态分布，但在实际的生产中，这个分布往往是未知的，本研究用总体观测值的最大熵分布近似估计其真实分布，基于此对马尔科夫链中的转移概率进行计算，并依次计算CUSUM控制图的ARL值，实现对CUSUM的评价。

3 仿真案例

为了详细地说明最大熵分布对Shewhart控制图构建及对CUSUM控制图评价方法的改进，本文设计了仿真案例来进行说明。在机械制造、电气等领域的质量控制中，考虑到质量特性参数服从重尾分布的情况比较多[33-34]，基于此，我们以威布尔分布为例，对控制图的性能进行分析。假设质量特性参数X的真实分布为威布尔分布，其概率密度函数f(X)为

(15)

其中λ为比例系数，k为形状系数，且λ>0，k>0。

由于在实际生产过程中，很难确切地知道质量特性参数的真实分布情况，本节首先应用最大熵估计对真实分布进行近似，然后依此构建最大熵控制图，并对控制图进行评价，最后基于最大熵分布应用马尔科夫链法计算自适应CUSUM控制图ARL值，并与真实分布下控制图的ARL值进行比较。

3.1 最大熵分布对威布尔分布的拟合

为了求解最大熵函数，我们利用Matlab软件进行仿真。生成1000个服从比例参数λ=1，形状参数k=2的威布尔分布的离散随机数，并把这些信息作为质量特性参数X最大熵分布的先验信息，参考Zhang和Xu[35]的研究，把X的先验信息转化为函数：u1(x)=lnx，u2(x)=x2因此根据公式(1)，可得到以下最大熵分布的规划模型：

(16)

从而得到最大熵分布函数为：

f(x)=exp[λ0+λ1ln(x)+λ2x2]

通过对约束条件中参数的计算，得到M1=-0.2589,M2=0.1968。利用Matlab软件的非线性规划工具，求解未知参数的值，得到：

λ0=0.6931，λ1=1.0022，λ2=-0.9874

因此，最大熵函数为：

f(x)=exp[0.6931+1.0022ln(x)

-0.9874x2]

(17)

为了验证最大熵函数的有效性，将求得的f(x)与威布尔分布的密度函数进行对比，如图3所示。从图中可以看出，此处最大熵分布对威布尔分布的拟合较好。因此基于最大熵分布构建Shewhart控制图和对CUSUM控制图进行评价，可以有效实现对质量特性参数的控制。

图3 最大熵函数与威布尔分布概率密度函数对比

3.2 基于最大熵分布控制图控制界限的确定

确定了质量特性参数X的最大熵分布函数f(x)之后，就可以确定Shewhart控制图和CUSUM控制图的控制界限。

基于最大熵分布对Shewhart控制图的改进。给定第Ⅰ类错误α=0.05，将公式(17)中的最大熵函数带入公式(9)中，通过数值法近似求得最佳控制区域为：A*=[0.076, 1.766]。而在传统的Shewhart控制图中，X被假设为服从正态分布，其控制图的控制界限设置为[μ-3σ,μ+3σ](通过对离散值的估计，得到均值和标准差分别为μ=0.8977，σ=0.4436)，得到X的Shewhart控制界限为[-0.433,2.228]。表1和图4给出了应用最大熵分布对Shewhart控制图进行改进前后的比较，从实际的威布尔分布看，X>0，而传统的控制图LCL的值小于0，因而没有意义。而改进后的控制图的控制区域比传统方法得到的控制区域更短，且更符合分布的实际情况。

图4 改进前后的Shewhart控制图的控制区域

表1 改进前后Shewhart控制图的控制界限

3.3 基于ARL的最大熵Shewhart控制图评价

用失控状态下的ARL来比较最大熵对Shewhart控制图进行改进前后的控制性能，即在第Ⅰ类错误α=0.05下，当过程失控时，ARL越小越好。而当过程失控时，X不服从公式(15)中的威布尔分布，因此本文用参数λ的偏移来表示过程失控，若参数λ的偏移数值越大，代表过程失控程度越严重。由于在X的原分布中，比例参数λ=1，形状参数k=2，则原分布概率密度函数为：

f(x)=2x·exp(-x2)

(18)

设过程失控时，即λ≠1。设失控状态下密度函数为f'(X)，则：

(19)

在参数λ不同的偏移下，通过f'(X)可计算第Ⅱ类错误β，即：

(20)

其中UCL和LCL分别为用改进前后控制界限计算方法得到控制区域的上控制限和下控制限(如表1所示)。

对于Shewhart控制图，根据公式(20)，可以求得第Ⅱ类错误β，并将其代入公式(11)，就可以得到改进前Shewhart控制图和改进后最大熵Shewhart控制图在失控状态下的ARL(如表2所示)。

表2 失控状态下改进前后Shewhart控制图的ARL值

由表2可知，随着观测值漂移量的增加，即λ的增大，最大熵Shewhart控制图和Shewhart控制图的ARL值都在减少。当偏移量较小时(λ<2，均值偏移量<2σ)，最大熵Shewhart控制图比传统Shewhart控制图的ARL小的多，因此在这种情况下，最大熵Shewhart控制图比传统Shewhart控制图能更快地检测到观测值的漂移。但当偏移量较大时(λ>2，均值偏移量>2σ)，改进后的Shewhart控制图ARL值也小于改进前的ARL值，因此也可以更快地发现过程偏移。总的来说，最大熵Shewhart控制图具有更好的监控效果。

3.4 基于最大熵分布的CUSUM控制图评价

为了验证基于最大熵分布的CUSUM控制图的评价方法。本文应用自适应CUSUM控制图对过程进行监测。为求控制界限，我们给出不同的偏移δ，取值为0.25～1.5，取值间隔为0.25。根据公式(12)，得到了不同偏移下，监测上偏移的自适应CUSUM控制图的控制界限，如表3所示。

表3 自适应CUSUM控制图的控制界限

对自适应CUSUM控制图的评价。本研究应用2.3.2中的马尔科夫链法计算控制图的ARL，在这个过程中，分别应用真实威布尔分布(公式(18))、最大熵分布(公式(17))以及正态分布(数学期望和标准差为：μ=0.8977，σ=0.4436)作为总体分布的估计，并进行转移概率的计算，最后得到三种分布下的ARL值，如表4所示。

表4 失控状态两种CUSUM控制图的ARL值

从表4可以看出，在小偏移的情况下(λ<2)，自适应CUSUM控制图的ARL要比相同偏移下的Shewhart控制图的小(对比表3)，这也说明了自适应CUSUM控制图对小偏移的监测效果比Shewhart控制图要好。而随着偏移增大，所有控制图的ARL都在减小。且当偏移后的比例参数λ为1.1至1.7时，基于正态分布 CUSUM控制图的ARL与真实值的差距较大，而基于最大熵分布计算的CUSUM控制图的ARL与真实值的差距较小。当λ超过2时，CUSUM控制图的ARL下降速度明显变慢。同时可以明显地看出，在威布尔分布下，基于最大熵的CUSUM控制图ARL值较正态分布下的结果更接近于真实值。这说明在重尾分布的情况下，应用基于最大熵的方法对CUSUM控制图进行评价的结果更接近与真实分布的情况，结果也更客观。

3.5 小结

从仿真案例的结果看，不同控制图在过程的不同偏移下监控效果不同，在实际情况中，应结合最大熵Shewhart控制图和自适应CUSUM控制图特点，即在过程发生不同偏移时采用不同的控制图进行监控。当过程发生的偏移较小时，可采用自适应CUSUM控制图进行监控；当过程发生的偏移较大且不一定服从正态分布时，应采用最大熵Shewhart控制图进行监控。当过程发生的偏移较大且服从正态分布时，可直接选择Shewhart控制图进行监控；在对CUSUM控制图的控制效果的评价中，基于质量特性参数最大熵分布的控制图评价方法能得到更接近于真实分布的结果。

4 结语

本文对质量特性参数为一般分布下的Shewhart控制图进行了改进和评价，并针对自适应CUSUM控制图，提出了基于最大熵的评价方法。首先，根据最大熵原理，对分布进行了拟合，求出了质量特性参数的最大熵分布，以此构建了最大熵Shewhart控制图，并选用失控状态下ARL作为控制图的评价指标；其次，针对自适应CUSUM控制图，本研究提出了结合最大熵的马尔科夫链方法对控制图性能进行评价。最后通过仿真实验，选取一个重尾分布进行验证。结果表明，基于最大熵Shewhart控制图的控制效果要优于传统的Shewhart控制图。基于最大熵的CUSUM控制图性能评价方法得到的结果与真实分布的情况更为接近。同时，本文还总结了不同偏移和分布情况下控制图的适用情况，为实际生产中进行质量控制提供参考。