1. 引言
随着科技的飞速发展,许多制造业、非制造业领域中的事件(如不合格品、缺陷)发生率低至百万分之一,这样被监测事件发生率极低的生产过程被称为高质量过程 [1] ,高质量过程中既定单位时间内事件的平均发生次数几乎为零,所以被观测事件的发生呈现低概率特性,这使得基于传统属性控制图(c控制图、p控制图、np控制图)存在高误报率、低检测效率等风险。幸运的是,事件的低发生概率恰恰对应大间隔时间,通过转换监测特征尺度,以两个相邻事件发生的间隔时间(Time Between Events, TBE)为质量特征,可以克服传统控制图缺陷,提高高质量过程监控性能 [2] 。另一方面,由于TBE的降低、增加分别意味着过程可能发生改进、恶化,而过程的恶化、改进与保证产品质量、改良生产设备等密切相关,所以及时监测过程的恶化、改进十分有必要。因此,开发基于TBE的双边控制图,以实现及时监测过程的改进、恶化情况 [3] 。
在事件发生率恒定的情况下,TBE可以用指数分布很好建模。而指数分布具有偏态性特点,这使得TBE不能直接应用基于正态分布传统控制图进行监测 [4] [5] 。Chan等 [6] [7] 先后提出了累积质量(Cumulative quantity control, CQC)控制图、累积概率控制(Cumulative probability control, CPC)图;C. W. Zhang等 [8] 发现这些基于传统概率极限的指数图存在平均运行长度(Average Run Length, ARL)偏差问题。Johnson [9] 指出若X服从尺度参数为的
指数分布,则变换后的变量
服从形状参数
和尺度参数
的威布尔分布。Nelson、Liu [10] [11] 经分析发现幂次为3.6或2.5的变换,能够使得指数分布下的随机变量近似正态分布。2005年Montgomery [12] 解释了利用变量变换构建TBE控制图的思想并建议利用在小偏移监测方面表现良好的累计和(Cumulative sum, CUSUM)、指数加权移动平均(Exponentially weighted moving average, EWMA)控制图来进行监测。McCool等 [13] 证明了恰当的幂变换可以有效地提高EWMA控制图的性能。Santiago和Smith (2013) [14] 提出了添加运行规则后的休哈特控制图,进一步,为监测到小偏移,Aslam等 [15] 提出了变量变换下双边的双指数移动平均(double exponentially weighted moving average, DEWMA)控制图,发现其控制图的性能优于文献 [13] 。考虑到正态分布下DEWMA控制图比 EWMA控制图更灵敏的优点,Adeoti [5] 提出了基于变量变换的双边DEWMA控制图来监测遵循指数分布的质量特性,以改进文献 [15] 提出的EWMA控制图。
在高质量过程中,尽早发现过程恶化或过程的改进对减少生产过程中带来的损失或降低生产成本具有重要的意义。另一方面,考虑到高质量过程中突发事件可能引发灾难性损失,而这类事件会导致样本的最新差异信息陡变。Naveed (2018)等 [16] 提出的扩展指数加权移动平均(Extended Exponentially Weighted Moving Average, EEWMA)控制图,其是EWMA控制图的一种扩展形式,其利用了当前、过去的信息和当前的最新变化信息,实验证明,相对于EWMA控制图,EEWMA控制图对正态分布下的小偏移更具有灵敏性。
此外,注意到文献 [5] 、 [15] 分别利用概率的方法计算了所提控制图的ARL。Maragavio、Adams等 [17] [18] 指出,采用概率方法计算ARL不适用于CUSUM、EWMA等记忆型控制图。Haq等 [19] 指出,针对EWMA控制图以及其一些改进形式,在假设信号事件可以由一系列独立的伯努利随机变量表示下近似计算的ARL可能非常不准确,并认为需要使用马尔可夫链或计算机仿真模拟来重新评估。
综上,本文针对高质量过程中TBE的监控问题,结合EEWMA控制图的优点,并以指数分布下TBE作为过程质量的量化指标。通过幂次为3.6的幂变换,实现数据近似正态化,开发了指数分布下基于幂变换的EEWMA控制图,简称为PT-based EEWMA控制图。考虑到实际生产中事件发生率可能变化,因此分析了所提控制图在威布尔分布下的稳健性。同时,通过蒙特卡洛模拟重新评估了指数分布下基于3.6幂变换的EWMA和DEWMA控制图的ARL性能,并与变量变换后的PT-based EEWMA控制图进行比较。
2. 指数分布下基于数据幂变换的EEWMA控制图的设计
2.1. 现有正态分布下的双边EEWMA控制图
在本节我们简要介绍Naveed [16] 提出的正态分布下监控均值偏移的EEWMA控制图。假设随机样本
,
独立且服从正态分布,过程受控时均值、方差分别为
、
。EEWMA控制图统计量定义为
, (1)
其中,
和
为平滑系数,
,
,设
,“
”表示当前最新的变化信息。整理(1)式可得
, (2)
EEWMA控制图统计量的均值和方差分别为
, (3)
, (4)
其中,
。当
时,EEWMA控制图简化为具有时变控制限的EWMA控制图。EEWMA控制图的控制限表示为
, (5)
, (6)
k为控制系数,其值可通过指定的受控平均运行长度ARL0确定。
2.2. 指数分布下基于幂变换的EEWMA控制图的设计
在高质量过程中,将事件间隔时间TBE记为随机变量X,则
,其概率密度函数为
, (7)
其中,尺度参数
,X的期望
和方差
分别为
和
。设
,
表示过程处于受控状态,
表示过程处于失控状态。
为X独立同分布的随机样本。
已知均值为
的指数分布是尺度参数为
、形状参数为1时的威布尔分布,即
。并且经过幂变换
后的变量Y仍服从威布尔分布,其形状参数和尺度参数分别为
和
,根据文献 [10] ,经过
幂变换后的数据能够近似正态分布,因此,令
, (8)
其均值和方差分别为
, (9)
, (10)
构建PT-based EEWMA控制图的统计量
(11)
,
,令
。PT-based EEWMA控制图统计量的均值和方差分别
, (12)
, (13)
基于式(12)、(13),构建在t点处的控制限
其中,
表示控制图系数,当
或
时,控制图在t点发出失控信号。
3. 性能评价
通常利用平均运行长度(ARL)对控制图的性能进行评价,ARL是控制图从开始监测到失控信号发出的TBE平均观测样本数,ARL0指过程实际处于受控状态时的平均运行长度,所以ARL0越大越好,使得控制图具有尽可能小的误报率;而ARL1是过程处于失控状态时的平均运行长度,我们希望其越小越好,以保证控制图尽早检测到过程的偏移。
Margavio [17] 证明了在假设运行长度RL服从参数为p的几何分布背景下,利用
计算ARL的方法不适用于记忆性控制图,Abbasi [19] 、Haq [20] 等人也指出这种利用概率近似求EWMA及其一些改进形式的控制图的ARL方法存在问题。因此我们用蒙特卡洛模拟方法重新计算文献 [5] 、 [15] 所提控制图的ARL,并与新提出的控制图进行比较分析。
3.1. 现有控制图
本小节主要对文献 [5] 、 [15] 提出的指数分布下基于3.6幂变换的EWMA、DEWMA控制图进行简单介绍。
3.1.1. 现有基于幂变换的EWMA控制图
Alsam等人 [15] 首先对指数分布的数据进行幂次为3.6的幂变换,基于式(8)经变换得到的
,定义在t时间点的图统计量
, (14)
其中平滑系数
,统计量的均值和方差分别为
,
,控制限为
, (15)
, (16)
其中,
为控制图系数,当
或
时,控制图发出失控信号。
3.1.2. 现有基于幂变换的DEWMA控制图
Adeoti等人 [5] 开发了对指数分布数据进行幂次为3.6的幂变换下的DEWMA控制图,基于式(8)经变换得到的
,图统计量定义为
, (17)
文献 [5] 分别对
、
下的控制图性能进行了比较分析,结果表明,
时控制图具有更好的监测性能。因此本文关注
时控制图的性能,设
,此时图统计量的均值和方差分别为:
,
控制限为:
其中,
为控制图系数,当
或
时,控制图发出失控信号。
3.2. 控制图的性能比较
本节以ARL作为衡量指标,利用Monte Carlo模拟10,000次,指定
,选取PT-based EEWMA控制图中
,对应
在
的范围内取不同的值,以实现在不同平滑系数下对控制图监测性能的全面分析,表1~4分别展示了新提出的PT-based EEWMA控制图与文献 [5] 、 [15] 提出的指数分布下基于3.6幂变换的EWMA、DEWMA控制图在不同平滑系数和不同大小、方向偏移下的ARL1值。

Table 1. When λ 1 = λ = 0.05 , the ARL of EEWMA, EWMA, and DEWMA control chart
表1.
时,数据幂变换下EEWMA、EWMA、DEWMA控制图的ARL值

Table 2. When λ 1 = λ = 0.1 , the ARL of EEWMA, EWMA, and DEWMA control chart
表2.
时,数据幂变换下EEWMA、EWMA、DEWMA控制图的ARL值

Table 3. When λ 1 = λ = 0.3 , the ARL of EEWMA, EWMA, and DEWMA control chart
表3.
时,数据幂变换下EEWMA、EWMA、DEWMA控制图的ARL值

Table 4. When λ 1 = λ = 0.5 , the ARL of EEWMA, EWMA, and DEWMA control chart
表4.
时,数据幂变换下EEWMA、EWMA、DEWMA控制图的ARL值
模拟结果显示(见表1~4):
1) 对于较大尺度的上、下偏移,
固定时,
越小,PT-based EEWMA控制图的ARL1越小。而对于较小尺度的上、下偏移,相对较大的
对应的控制图ARL1较低。此外,当
固定时,
越小,PT-based EEWMA控制图的ARL1越小。例如在k = 0.5的偏移下,
固定为0.03,
取0.05、0.1、0.3、0.5对应的ARL1值分别为:18.11、21.32、50.56、107.46。
2) 对比指数分布下基于幂变换的EWMA控制图,除了当
对应的PT-based EEWMA控制图在个别较大偏移下(如k = 0.1、5、7)的ARL1略高于该EWMA控制图以外,对于其他不同大小的上下偏移,不同的平滑系数下的PT-based EEWMA控制图ARL1均低于该EWMA控制图。
3) 对比指数分布下基于幂变换的DEWMA控制图,当
取值较小时,PT-based EEWMA控制图在大部分偏移下的ARL1值都低于该DEWMA控制图(如
),对于较大的
取值(如
),PT-based EEWMA控制图在TBE均值的下偏移监测方面不如该DEWMA控制图表现好,但在TBE均值发生上偏移时,PT-based EEWMA控制图的性能均优于该DEWMA控制图。总的来看,较小的
值对应的PT-based EEWMA控制图的性能优于指数分布下基于幂变换的DEWMA控制图,并且当
取值较大时,在监测下偏移方面,控制图也具有更好的监测性能。
总的来说,针对TBE均值不同大小、方向下的偏移,PT-based EEWMA控制图整体上表现更好,并且两个平滑系数的引入使得控制图能够更加灵活地用于对不同大小偏移的监测。
4. PT-based EEWMA控制图对威布尔分布的稳健性分析
上述内容是基于事件发生率不变,即事件间隔(TBE)服从指数分布的假设下进行的研究,鉴于实际高质量生产过程中,事件发生率可能随着时间发生变化,此时威布尔分布是拟合TBE更好的模型。Liu等 [11] 分析了当TBE真实分布为威布尔分布时,指数分布下基于幂变换的EWMA图的稳健性。在本节我们假设事件之间的时间服从指数分布,分析当真实分布为威布尔分布时,指数分布下PT-based EEWMA控制图的稳健性。
威布尔分布
的密度函数如下所示:
其中,
,
是尺度参数,
为形状参数。当
时,威布尔分布简化为指数分布。设
,已知经过
变换后,Y仍服从威布尔分布,即
,其均值和方差分别为:
文献 [11] 利用马尔可夫链方法分析了对指数分布数据进行双平方根变换后的EWMA控制图的ARL性能,发现其受控ARL仅与平滑系数、控制图系数有关,与尺度参数的大小无关,且此结论对真实分布为威布尔分布的情况也成立,即当形状参数固定不变时,尺度参数的变化不影响其受控ARL的值。
在本文我们利用蒙特卡洛方法验证了对指数分布TBE数据进行幂次为3.6的幂变换后的PT-based EEWMA控制图来说,上述性质同样成立。PT-based EEWMA控制图的受控ARL值仅与平滑系数
、
以及控制图系数的取值有关,不受指数分布下尺度参数变化的影响,且当真实分布为威布尔分布且形状参数不变时,控制图的受控ARL不受尺度参数变化的影响。然而在威布尔分布下,形状参数的变化会导致受控ARL发生变化。不失一般性,假设
,展示当真实分布为威布尔分布时,不同的
和
平滑参数下形状参数
变化对控制图的受控ARL性能的影响,详细见表5
从表5中可以发现,
越小,控制图的稳健性越好。并且当
固定时,
越小,控制图的稳健性越好,如
时,对于(0.6, 4)范围内的
,受控ARL始终在500的3%以内。

Table 5. When the real distribution is the Weibull distribution, the ARL0 of the EEWMA control chart based on the data power transformation
表5. 真实分布为威布尔分布时,数据幂变换下EEWMA控制图的ARL0
表6、表7分析了在不同形状参数以及不同平滑系数选取下,面对不同大小、方向的偏移时控制图的性能。结果表明,形状参数
的取值对失控情况下控制图的监测性能具有显著影响。对于不同大小的偏移,
时控制图对偏移的敏感度优于均比
,此外,
时存在极大ARL1值的现象,如
,
,
时,ARL1的值达到了1471.88。此外,随着形状参数的增大,控制图的监测性能越好。

Table 6. When the real distribution is the Weibull distribution, the ARL1 of the EEWMA control chart based on the data power transformation ( λ 1 = 0.1 )
表6. 真实分布为威布尔分布时,数据幂变换下EEWMA控制图的ARL1 (
)

Table 7. When the real distribution is the Weibull distribution, the ARL1 of the EEWMA control chart based on the data power transformation ( λ 1 = 0.3 )
表7. 真实分布为威布尔分布时,数据幂变换下EEWMA控制图的ARL1 (
)
5. 应用
在本节中,我们运用模拟数据集以及与尿路感染(UTI)相关的真实数据集来展示PT-based EEWMA控制图的实施过程,验证了PT-based EEWMA控制图相比于现有指数分布下基于幂次为3.6幂变换的EWMA控制图具有更好的监测性能。
5.1. 实例应用
利用构建的PT-based EEWMA控制图来实施对尿路感染(UTI)情况的监测,数据是来自文献 [14] 使用的来自大型医院系统的真实数据,此数据也先后被文献 [5] [15] 使用,该医院着重关注尿路感染的高发率,并通过监测在医院期间尿路感染患者的出院频率,以及时发现感染率是否有上升或下降的趋势。本节针对男性进行分析,数据见表8,已知男性UTI患者之间的平均时间为0.21天,我们可以假设受控状态的
。不失一般性,指定
,选取PT-based EEWMA控制图平滑系数
,
,在选取指数分布下基于幂变换的EWMA控制图的平滑系数
,其控制图系数分别为2.688、2.687下,讨论我们的方法并与EWMA控制图进行比较分析,详细见图1。

Table 8. Urinary tract infection data
表8. 尿路感染(UTI)数据

Figure 1. EWMA and EEWMA control charts based on the UTI data power transformation for monitoring data in Table 8
图1. 幂变换下用于监控表8中UTI数据的EWMA、EEWMA控制图
图1展示了利用PT-based EEWMA控制图、基于幂变换的EWMA控制图对UIT进行监测的实施过程,我们可以看到结果与原文献一致,即处于受控状态,可将PT-based EEWMA控制图用于对将来过程的监控。
5.2. 模拟应用
本节通过模拟数据来分析控制图的失控性能,假设两个连续事件之间的时间服从
的指数分布,利用R软件生成54个来自
指数分布下的模拟数据,将其视为来自失控过程的TBE观测值,即过程中发生了向上的偏移。实施基于幂次为3.6的幂变换下的EWMA、PT-based EEWMA方案,指定
,选取PT-based EEWMA、EWMA控制图平滑系数分别为
,
,详细的实施过程见图2

Figure 2. EWMA and EEWMA control charts based on the data power transformation for monitoring simulated data
图2. 幂变换下用于监控模拟数据的EWMA、EEWMA控制图
如图2所示,可以发现PT-based EEWMA控制图与幂变换下的EWMA控制图在监测性能方面存在差异,PT-based EEWMA控制图在第15个点处的值超出控制限,即发出失控信号,幂变换下的EWMA控制图在第21个点处发出失控信号,说明PT-based EEWMA控制图更早监测出偏移,即PT-based EEWMA控制图的监测性能优于幂变换下的EWMA控制图。
6. 结论
在高质量过程中,尽早发现过程的恶化或改进对减少生产过程中带来的损失或降低生产成本具有重要的意义。本文对服从指数分布的TBE进行幂次为3.6的幂变换使其近似服从正态分布,结合EEWMA控制图引入两个平滑系数、利用最新差异信息使控制图更加灵活且对中、小偏移更加敏感的优点提出了针对高质量过程监控的PT-based EEWMA控制图。数值实验验证了相比于EWMA,PT-based EEWMA控制图对各种尺度的均值偏移均具有较优的ARL性能。此外,在
值较小时,PT-based EEWMA控制图对上下偏移的敏感度优于相应的DEWMA控制图,并且当
取值较大时,PT-based EEWMA控制图在监测下偏移方面也表现出更好的监测性能。最后进一步研究了事件发生率随时间变化因而TBE服从威布尔分布的情形。数值实验表明,具有较小平滑系数的PT-based EEWMA控制图,对于服从威布尔分布的TBE数据亦具有相当的稳健性。
基金项目
长安大学中央高校基本科研业务费专项资金资助项目(310812163504)。
NOTES
*通讯作者。