一种面向三元空间相关Poisson过程的控制图设计

吴苍, 司书宾

1.西北工业大学 机电学院, 陕西 西安 710072; 2.兰州理工大学 机电学院, 甘肃 兰州 730050

生产过程中对异常因素的识别和诊断是维持过程稳定、安全运行的重要环节。过程监控方法目前已从生产制造领域发展到服务、医疗和交通等领域,对过程中的异常防微杜渐,能够有效避免大的质量和安全问题。近年,随着传感技术和计算机技术的发展,对复杂过程安全性、可靠性和高效性能的监控提出了更高的要求。例如,在缺陷检测方面,具有多种质量特性的工业产品中经常出现不同类型的缺陷,这些缺陷之间有可能相互作用[1];在医疗领域,某一时间段内,相邻区域感染某种疾病的人数存在相关性[2]。

控制图是统计过程控制领域的一种重要工具。通过对过程数据的采集、分析、计算和打点,实时监控过程的运行状况,并对过程异常及时报警,从而达到监控的目的。根据过程数据的维度,其分为单值控制图和多元控制图。常见的单值控制图有适用于监控连续变量的休哈特图(Shewhart chart)、EWMA (exponentially weighted moving average)图和CUSUM (cumulative sum)图;监控一元离散变量的有c图、u图等。这些控制图应用范围很广,但是其假设观测值是独立同分布的,直接将其运用到多元情形会忽略数据间的关联性。为此，多元控制图应运而生,主要包括用于监控离散变量的np图、D图[3]，用于监控连续变量的多元Shewhart图(T2图)、多元EWMA图[4]和多元CUSUM图。其中,np图和D图将多元离散数据转化成一元数据,忽视了多元数据的相关性;多元Shewhart图(T2图)、多元EWMA图和多元CUSUM图不但研究了多个变量,并且同时考虑了变量之间的相关性,但这些控制图假设变量服从正态分布。综上所述,这些控制图在多元离散过程监控中有局限性。

因此有必要建立一种针对多元离散型过程数据的控制图。最初,Niaki和Abbasi[5]考虑把多元Poisson分布转化为近似的多元正态分布进而采用T2控制图监控。但当数据的均值较小或者存在很多零值时,这种近似会产生较大误差。Chen等[6]假设Poisson过程的参数服从多元对数正态分布,并构建多元EWMA控制图监控多元Poisson过程,但这一模型忽视了不同维度数据间的空间关联性,在应用中有一定的局限性。He和龙威等[7-8]分别在等方差和异方差多元Poisson模型基础上,构建了CUSUM控制图。这类模型的相关性度量不够灵活,不适应如非线性或对称性的关联结构。

由于copula模型在处理非线性、非正态等结构上的优势,在金融和经济领域得到广泛应用。与此同时,这一方法也被引入到多元数据建模和监控领域。Dokouhaki等[9]基于copula函数构建面向二元自相关二项型数据的CUSUM控制图。Sasiwan-napong等[10]假设观测值服从指数分布,应用多种copula构建T2和多元EWMA控制图。Sukparungsee等[11]采用5种copula(Normal,Clayton,Frank,Gum-bel,Joe copula)构建T2控制图用于监控指数型特征值。Kuvattana等[12]针对二元指数型数据构建了多元EWMA和多元CUSUM控制图,仿真结果显示多元EWMA控制图的监控效果最好,并且基于Normal copula构建的多元EWMA图可以用于监控正负相关及各种强度相关性下的漂移。针对多零值离散数据,Fatahi等[13]应用零膨胀Poisson分布(zero inflated Poisson distributions)构建基于copula的联合概率分布监控二元稀少事件(rare event)。考虑到非正态变量,Krupskii等[14]构建了对密度函数形状变化敏感的多元过程监测方案。Chen[15]采用Normal copula描述时间序列的自相关和多元数据的互相关,并实施了多元EWMA控制图。近期,Song等[16]提出基于copula的非参数方案,监控二元未知分布的位置和规模信息。综上所述,多数研究聚焦于二元或连续型数据建模及监控。这主要有两方面原因：①把二元copula函数推广到多元情形会产生较复杂的结构,也给参数估计带来困难；②离散边缘分布往往会产生不唯一的copula函数。因此,有必要探索针对三元离散数据的建模和控制图研究。

本文以三元Poisson过程数据为对象,考虑各时刻不同维度数据的空间相关性,引入三元Clayton copula描述各维度数据的相关性,根据对数似然比检验理论,设计出与模型匹配的控制图。采用马尔科夫链法近似求解出多种参数变化下的平均运行链长,最终验证控制图的监控性能。

1 Clayton copula简介

1.1 copula函数简介

copula模型广泛应用于描述和度量变量间的非线性的相关结构[17]。以常见的二元copula函数为例,二元变量(y1,y2)的边缘分布函数和联合分布函数为F1,F2和F。简便起见,记v1=F1(y1),v2=F2(y2)，且v1∈[0,1],v2∈[0,1],copula函数记做C(v1,v2),存在如(1)式所示关系

F(y1,y2)=C(F1(y1),F2(y2))=C(v1,v2)

(1)

式中，C(v1,v2)∈[0,1],且满足以下性质:

1) 对于所有的v1,v2∈[0,1],存在

如果边缘分布函数F1,F2均连续,那么copula函数是唯一的,反之,copula函数往往不唯一。限于此,copula函数目前在连续型过程得到充分应用,而在离散和混合过程应用较少。

根据构造原则常常把copula分为阿基米德族和椭圆族。阿基米德函数族copula结构类似,根据生成函数的不同将其分为Clayton,Frank和Gumbel函数等。而椭圆族copula是由椭圆族函数推导而来,常见的有Gaussian和Student-t等。又根据参数的数量分为单一参数copula和多参数copula。各种类型的copula在描述相关性时有所差别[18]。Gaussian copula可以灵活地描述同等程度的正相关和负相关,并在其允许范围内包含Fréchet边界[18]。Clayton copula不适用于解释负相关关系,但适用于高强度的正向相关。Frank copula在实践中得到广泛应用,因为Fréchet边界包含在其允许范围内。

1.2 Clayton copula模型

Clayton copula是阿基米德函数族中常见的单一参数模型[17],广泛应用于金融和经济领域。二元Clayton copula定义为

(2)

式中，β∈(-1,+∞){0}。Clayton copula是常见的单变量copula模型,生成函数为

满足φ(C(v1,v2))=φ(v1)+φ(v2)性质。

三元Clayton copula[19]定义为

(3)

式中,β>0是唯一参数,表示任意2种变量(v1和v2;v1和v3;v2和v3)的相关性。该模型有很强的对称性并且结构简单,适用于三元数据建模。

2 基于copula的三元Poisson模型

实际生产过程中,多种质量特性在同一时刻的采样值往往存在相关信息。某一时刻t采集到M个离散的观测值,记做Yt=(Yt,1,…,Yt,M)。假设Yt服从多元Poisson过程,不同时刻采集的数据间相互独立(如Yt,Yt-1),且每个时刻的观测值存在空间相关性(如Yt,1,Yt,2)。考虑Yt的空间相关性,用copula来构建其联合概率质量函数。假设对所有t=1,…,T,向量Yt=(Yt,1,…,Yt,M)的联合分布为

F(Yt;Θ)=C(F1(Yt,1),F2(Yt,2),…,FM(Yt,M);β)

(4)

式中：β代表copula参数；Fi代表边缘累积分布函数,i=1,…,M。

对于二元Poisson过程,在时刻t获得的2个观测值可以表达为向量Yt=(Yt,1,Yt,2),其联合分布和概率质量函数分别表示为

(5)

同理,三元Poisson过程中变量Yt=(Yt,1,Yt,2,Yt,3)的联合分布和概率质量函数分别表示为

(6)

三元Poisson过程的协方差矩阵表示为

(7)

式中，apq=aqp且apq=cov(Y*,p,Y*,q),p,q=1,2,3。根据协方差的定义结合copula理论可以计算

式中，E(Y*,p)是Y*,p边缘分布的均值,(Y*,p,Y*,q)的二元联合概率质量函数为

hpq(Y*,p,Y*,q;β,λp,λq)=

显然,对于Poisson分布每个维度变量的取值范围是非负整数集。由于各变量的协方差不恒定,本模型也属于异方差模型。

实践中,分布参数λ及copula参数β往往是未知的,需要根据一系列稳态观测值估计这些参数。两阶段最大似然估计法,也称作IFM(inference functions for margins)方法[20],是一种常用的多元copula模型的参数估计方法,表述如下:

(8)

(9)

两阶段估计法在连续情形效果显著,并得到了广泛应用。而在离散情形下,估计参数在小维度情形很有效,对于大维度数据的效果次之[20]。

3 三元CUSUM控制图设计

针对多变量监控问题,联合监控多个变量比单独监控一个变量更加容易发现过程异常[19]。考虑到空间相关性建模并获得适用的受控模型及参数后,需要结合数据特点及控制图原理设计出适用的控制图。本文结合第2节的三元相关Poisson数据模型,设计多元CUSUM控制图。假设同一时刻的三元Poisson数据存在空间相关,而不同时刻的数据间相互独立。记t时刻的三元Poisson数据为Yt=(Yt,1,Yt,2,Yt,3),t=1,…,T,各个维度的参数记做λ=(λ1,λ2,λ3)。假设受控状态的Poisson参数为λ0，偏移发生后变化为λc，变点模型可表示为

根据似然比检验的思想,对变点存在和不存在2种情形的联合密度函数做似然比,构造统计量

(10)

式中,f(Ym,…,Yn;λ0)和f(Ym,…,Yn;λc)是(Ym,…,Yn)分别在λ=λ0和λ=λc时的联合密度函数。记

(11)

式中，f(Yi)是概率质量函数。对于本文所讨论的二元Poisson过程,f(Yi)=hB(yi,1,yi,2);同理,三元过程中,f(Yi)=hT(yi,1,yi,2,yi,3)。基于此,似然比统计量可以进一步写做

(12)

等价的递归形式可表示为

Sn=max{0,Sn-1+ln(rn)}

(13)

式中，S0=0。一旦Sn超出给定的控制限h则会报警,意味着存在变点m,即Poisson参数λ发生了改变;否则,认为过程参数λ没有显著改变,即过程处于受控状态。

4 仿真实验与分析

4.1 多元相关Poisson数据的生成

基于copula模型生成多元数据是建立在条件分布基础上的。由于copula模型的多样性,数据生成存在一定差异。本节以三元Poisson过程为例,生成基于Clayton copula的随机数据。记变量v1和(v1,v2)的条件概率密度函数分别为C1(v1)和C1(v1,v2),存在

另记Ck(vk;v1,…,vk-1)为条件分布,即给定(v1,…,vk-1)条件下vk的条件分布。对三元变量问题,k=1,2,3。若分子分母同时存在且分母不为零时

Ck(vk;v1,…,vk-1)=

(14)

当k=1时显然存在C1(v1)=v1,当k=2,3时可以推算出

生成一组三元Poisson数据的过程如下所示：

1) 从均匀分布U(0,1)中随机生成v1;

2) 从均匀分布U(0,1)生成随机数作为C2(v2;v1)的实际值,结合v1求解方程得v2;

3) 从均匀分布U(0,1)生成随机数作为C3(v3;v1,v2)的实际值,结合v1和v2求解方程得v3;

4) 应用Poisson分布的反函数得到一组数据(y1,y2,y3),其中,y1=F-1(vi),i=1,2,3。

4.2 参数估计

4.3 马尔科夫链近似求解平均运行链长

目前,主流的用于评价控制图监控效果的指标是平均运行链长(average run length,ARL),即运行链长(RL)的期望值。而运行链长RL是指从开始到第一次失控报警的样本点数目。本文采用马尔科夫链(Markov chain)近似计算多元CUSUM控制图平均运行链长。控制限记做h,受控的多元Poisson参数为λ0=(λ1,λ2,λ3),目标参数记做λc=θλ0,失控参数λd=(δ1λ1,δ2λ2,δ3λ3)。首先将受控区间分为ne等分,每一段的长度Δ=h/ne,每部分区间可表示为[(k-1)Δ,kΔ),k=1,2,…,ne。一旦统计量落入区间k,认为统计量近似取区间中点值,即(k-0.5)Δ。假定转移矩阵P中的每一个元素pjk代表统计量从区间j转移到区间k的概率,记做

pjk=P{Sn∈[(k-1)Δ,kΔ)|Sn-1=(j-0.5)Δ}

(17)

式中，j,k=1,2,…,ne。转移概率的计算首先需要列举出所有可能的三元Poisson数据组合,并分别计算其联合概率。需要注意的是,计算受控ARL(L0)时,采用hT(Yt;λ0,β);反之,计算失控ARL(L1)时,采用hT(Yt;λd,β)计算。对转移矩阵做一些变换即可得到表示ARL的向量L

L=(I-P)-11

(18)

式中,I是单位矩阵,1是全为1的列向量。由于初始值S0=0,因此向量L的第一个元素即为ARL的近似值。值得注意的是,该方法近似计算ARL的精确程度受到间隔的等分数ne影响。ne越大,计算的ARL越接近真实值,但是计算时间会相应增加。

4.4 三元CUSUM控制图的监控效果

在缺陷监测和疾病预防的背景下,本文关注多元Poisson数据的异常增长,采用单侧向上的控制图考察监控效果。同理,单侧向下或双侧控制图可以根据实际情况选择。鉴于多元Poisson数据的特点,控制图的控制下限不做要求,控制上限记做h。在阶段一,给定适当的受控ARL, 采用马尔科夫链的方法搜索控制限h。

对三元Poisson过程的每个维度分别考虑3种异常参数,则会有33=27种组合。显然,实现每一种组合的异常参数集合是不现实的。因此,考虑了9种组合方式。每一种组合方式中,参数变化幅度可记为δ=(δ1,δ2,δ3),又因为受控参数λ0=(λ1,λ2,λ3),则异常参数λd=(δ1λ1,δ2λ2,δ3λ3)。此外,控制图设计中的目标参数λc=(θλ1,θλ2,θλ3)。考虑到设计过程的便利性,设定为1.1,1.5和2.0。

表1 L0=200时多元CUSUM控制图在3种copula参数和3种目标参数下的控制限

设定λ0=(10,10,10)且L0=200,表1和表2分别展示了3种copula参数和3种目标参数条件下的控制限和失控平均运行链长。其中,序号1～3仅有1个维度的参数发生向上阶跃,4～6有2个维度发生变化,7～9表示3个维度均发生变化的情形。文献[4,7]采用一种表征马氏距离的变化尺度(shift size)量,或称为非中心性参数(non-centrality para-meter)

Ss=[(λc-λ0)Σ-1(λc-λ0)′]1/2

式中,λ0和λc分别是受控和失控的参数向量。如公式(7)所示,Σ是协方差矩阵。非中心性参数能够度量过程参数的整体变化。对于2组异常参数,如果其非中心性参数相等,则L1也相当。例如,表2中序号1代表的异常参数的变化幅度为δ=(1.1,1,1),其性能与δ=(1,1.1,1)及δ=(1,1,1.1)情形相当。 同理,序号6代表的异常参数的变化幅度为δ=(2,1,2),其性能与δ=(2,2,1)及δ=(1,2,2)情形相当。因此,表2所列出的9种组合方式足够代表全部27种异常情形的ARL性能。从表2中可以看出,当参数有变化时,L1值均小于L0值,即L0=200,说明本文方法可以检测多元Poisson过程均值向上偏移。并且,随着参数变化的幅度值不断增大,L1越来越小。对于一种特定的失控状态,随着copula参数的增加,平均运行链长普遍增大。也就是说,多元过程的空间相关性越来越强时,监控异常变化的难度也越来越大。对比3种目标幅度值来看,θ=1.1的控制图对δ1=δ2=δ3=1.1(序号7)的异常情况效果好;θ=1.5的控制图对δ1=δ2=δ3=1.5(序号8)效果好;θ=2的控制图对δ1=δ2=δ3=2(序号9)效果好。因此,设计多元Poisson控制图时,对较小的变化幅度需要选择较小的θ, 反之,对于可能存在的较大变化,选择较大的θ可获得更理想的效果。

表2 多元CUSUM控制图在3种copula参数和3种目标参数下的失控平均运行链长

图1展示了Poisson数据及其在多元CUSUM控制图的应用。设定异常参数的变化幅度δ=(1.1,1.1,1.1),则目标参数的变化幅度选为θ=δ。三元Poisson数据由4.1节方法生成。图1a)～1c)分别代表β=1,3,5时对应的多元Poisson数据,而对应的CUSUM统计量绘制在图1d)～1f)中。每组Poisson数据包括3个维度,每个维度有20个受控数据和20个失控数据,图中用蓝色虚线分隔。红色水平线代表在L0=200水平下的控制限。从表1中可得出在β分别取1,3,5时对应的控制限为2.5, 2.49和2.55。从图1中可以看出,对于所列出的3种β取值,提出的方法对失控数据都有较好的检出力。

图1 Poisson过程数据及其在CUSUM控制图的应用(L0=200)

图2对比了本文所提方法与文献[3]中D图的性能。为了便于对比,实验假定δ=δ1=δ2=δ3,其中δ=1代表受控状态,δ=1.1,1.2,…,2代表失控状态。图2a)～2c)分别展示在β=1,3,5条件下D图和本文所提多元CUSUM控制图在θ=1.1,1.5的ARL对比。从图中可以看出,D图在β=1条件下监控较大偏移量(θ≥1.2)性能优于本文所提CUSUM控制图。当β=3,5时,本文提出的CUSUM控制图对所有偏移量均取得更小的L1。其中θ=1.1的CUSUM控制图对较小的偏移量监测效果最好,θ=1.5的CUSUM控制图对较大的偏移量监测效果最好。这一点与表2的结论一致。

图2 D图与本文提出CUSUM控制图的性能对比

5 结 论

针对具有空间对称相关结构的三元Poisson过程监控问题,首先采用三元Clayton copula函数描述其相关性,并构建联合分布;在假设copula参数不变的前提下,基于对数似然比检验方法结合联合分布信息,推导出多元CUSUM控制图的监控统计量;在给定初始值的前提下,利用马尔科夫链方法近似求解ARL,验证其有效性,并与D图做对比。仿真结果表明,提出的方法能够有效监测过程均值的向上漂移,并且目标偏移量与实际偏移量相等时可以获得更好的监控性能。与D图相比,当数据间的相关性较强时,提出的方法对于所有的偏移量能取得更好监控效果。在未来工作中,将进一步研究时间和空间关联同时存在的多元数据建模及监控问题。