一種面向三元空間相關(guān)Poisson過程的控制圖設(shè)計

吳蒼, 司書賓

1.西北工業(yè)大學(xué) 機電學(xué)院, 陜西 西安 710072; 2.蘭州理工大學(xué) 機電學(xué)院, 甘肅 蘭州 730050

生產(chǎn)過程中對異常因素的識別和診斷是維持過程穩(wěn)定、安全運行的重要環(huán)節(jié)。過程監(jiān)控方法目前已從生產(chǎn)制造領(lǐng)域發(fā)展到服務(wù)、醫(yī)療和交通等領(lǐng)域,對過程中的異常防微杜漸,能夠有效避免大的質(zhì)量和安全問題。近年,隨著傳感技術(shù)和計算機技術(shù)的發(fā)展,對復(fù)雜過程安全性、可靠性和高效性能的監(jiān)控提出了更高的要求。例如,在缺陷檢測方面,具有多種質(zhì)量特性的工業(yè)產(chǎn)品中經(jīng)常出現(xiàn)不同類型的缺陷,這些缺陷之間有可能相互作用[1];在醫(yī)療領(lǐng)域,某一時間段內(nèi),相鄰區(qū)域感染某種疾病的人數(shù)存在相關(guān)性[2]。

控制圖是統(tǒng)計過程控制領(lǐng)域的一種重要工具。通過對過程數(shù)據(jù)的采集、分析、計算和打點,實時監(jiān)控過程的運行狀況,并對過程異常及時報警,從而達到監(jiān)控的目的。根據(jù)過程數(shù)據(jù)的維度,其分為單值控制圖和多元控制圖。常見的單值控制圖有適用于監(jiān)控連續(xù)變量的休哈特圖(Shewhart chart)、EWMA (exponentially weighted moving average)圖和CUSUM (cumulative sum)圖;監(jiān)控一元離散變量的有c圖、u圖等。這些控制圖應(yīng)用范圍很廣,但是其假設(shè)觀測值是獨立同分布的,直接將其運用到多元情形會忽略數(shù)據(jù)間的關(guān)聯(lián)性。為此，多元控制圖應(yīng)運而生,主要包括用于監(jiān)控離散變量的np圖、D圖[3]，用于監(jiān)控連續(xù)變量的多元Shewhart圖(T2圖)、多元EWMA圖[4]和多元CUSUM圖。其中,np圖和D圖將多元離散數(shù)據(jù)轉(zhuǎn)化成一元數(shù)據(jù),忽視了多元數(shù)據(jù)的相關(guān)性;多元Shewhart圖(T2圖)、多元EWMA圖和多元CUSUM圖不但研究了多個變量,并且同時考慮了變量之間的相關(guān)性,但這些控制圖假設(shè)變量服從正態(tài)分布。綜上所述,這些控制圖在多元離散過程監(jiān)控中有局限性。

因此有必要建立一種針對多元離散型過程數(shù)據(jù)的控制圖。最初,Niaki和Abbasi[5]考慮把多元Poisson分布轉(zhuǎn)化為近似的多元正態(tài)分布進而采用T2控制圖監(jiān)控。但當數(shù)據(jù)的均值較小或者存在很多零值時,這種近似會產(chǎn)生較大誤差。Chen等[6]假設(shè)Poisson過程的參數(shù)服從多元對數(shù)正態(tài)分布,并構(gòu)建多元EWMA控制圖監(jiān)控多元Poisson過程,但這一模型忽視了不同維度數(shù)據(jù)間的空間關(guān)聯(lián)性,在應(yīng)用中有一定的局限性。He和龍威等[7-8]分別在等方差和異方差多元Poisson模型基礎(chǔ)上,構(gòu)建了CUSUM控制圖。這類模型的相關(guān)性度量不夠靈活,不適應(yīng)如非線性或?qū)ΨQ性的關(guān)聯(lián)結(jié)構(gòu)。

由于copula模型在處理非線性、非正態(tài)等結(jié)構(gòu)上的優(yōu)勢,在金融和經(jīng)濟領(lǐng)域得到廣泛應(yīng)用。與此同時,這一方法也被引入到多元數(shù)據(jù)建模和監(jiān)控領(lǐng)域。Dokouhaki等[9]基于copula函數(shù)構(gòu)建面向二元自相關(guān)二項型數(shù)據(jù)的CUSUM控制圖。Sasiwan-napong等[10]假設(shè)觀測值服從指數(shù)分布,應(yīng)用多種copula構(gòu)建T2和多元EWMA控制圖。Sukparungsee等[11]采用5種copula(Normal,Clayton,Frank,Gum-bel,Joe copula)構(gòu)建T2控制圖用于監(jiān)控指數(shù)型特征值。Kuvattana等[12]針對二元指數(shù)型數(shù)據(jù)構(gòu)建了多元EWMA和多元CUSUM控制圖,仿真結(jié)果顯示多元EWMA控制圖的監(jiān)控效果最好,并且基于Normal copula構(gòu)建的多元EWMA圖可以用于監(jiān)控正負相關(guān)及各種強度相關(guān)性下的漂移。針對多零值離散數(shù)據(jù),Fatahi等[13]應(yīng)用零膨脹Poisson分布(zero inflated Poisson distributions)構(gòu)建基于copula的聯(lián)合概率分布監(jiān)控二元稀少事件(rare event)?？紤]到非正態(tài)變量,Krupskii等[14]構(gòu)建了對密度函數(shù)形狀變化敏感的多元過程監(jiān)測方案。Chen[15]采用Normal copula描述時間序列的自相關(guān)和多元數(shù)據(jù)的互相關(guān),并實施了多元EWMA控制圖。近期,Song等[16]提出基于copula的非參數(shù)方案,監(jiān)控二元未知分布的位置和規(guī)模信息。綜上所述,多數(shù)研究聚焦于二元或連續(xù)型數(shù)據(jù)建模及監(jiān)控。這主要有兩方面原因：①把二元copula函數(shù)推廣到多元情形會產(chǎn)生較復(fù)雜的結(jié)構(gòu),也給參數(shù)估計帶來困難；②離散邊緣分布往往會產(chǎn)生不唯一的copula函數(shù)。因此,有必要探索針對三元離散數(shù)據(jù)的建模和控制圖研究。

本文以三元Poisson過程數(shù)據(jù)為對象,考慮各時刻不同維度數(shù)據(jù)的空間相關(guān)性,引入三元Clayton copula描述各維度數(shù)據(jù)的相關(guān)性,根據(jù)對數(shù)似然比檢驗理論,設(shè)計出與模型匹配的控制圖。采用馬爾科夫鏈法近似求解出多種參數(shù)變化下的平均運行鏈長,最終驗證控制圖的監(jiān)控性能。

1 Clayton copula簡介

1.1 copula函數(shù)簡介

copula模型廣泛應(yīng)用于描述和度量變量間的非線性的相關(guān)結(jié)構(gòu)[17]。以常見的二元copula函數(shù)為例,二元變量(y1,y2)的邊緣分布函數(shù)和聯(lián)合分布函數(shù)為F1,F2和F。簡便起見,記v1=F1(y1),v2=F2(y2)，且v1∈[0,1],v2∈[0,1],copula函數(shù)記做C(v1,v2),存在如(1)式所示關(guān)系

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

(1)

式中，C(v1,v2)∈[0,1],且滿足以下性質(zhì):

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

如果邊緣分布函數(shù)F1,F2均連續(xù),那么copula函數(shù)是唯一的,反之,copula函數(shù)往往不唯一。限于此,copula函數(shù)目前在連續(xù)型過程得到充分應(yīng)用,而在離散和混合過程應(yīng)用較少。

根據(jù)構(gòu)造原則常常把copula分為阿基米德族和橢圓族。阿基米德函數(shù)族copula結(jié)構(gòu)類似,根據(jù)生成函數(shù)的不同將其分為Clayton,Frank和Gumbel函數(shù)等。而橢圓族copula是由橢圓族函數(shù)推導(dǎo)而來,常見的有Gaussian和Student-t等。又根據(jù)參數(shù)的數(shù)量分為單一參數(shù)copula和多參數(shù)copula。各種類型的copula在描述相關(guān)性時有所差別[18]。Gaussian copula可以靈活地描述同等程度的正相關(guān)和負相關(guān),并在其允許范圍內(nèi)包含F(xiàn)réchet邊界[18]。Clayton copula不適用于解釋負相關(guān)關(guān)系,但適用于高強度的正向相關(guān)。Frank copula在實踐中得到廣泛應(yīng)用,因為Fréchet邊界包含在其允許范圍內(nèi)。

1.2 Clayton copula模型

Clayton copula是阿基米德函數(shù)族中常見的單一參數(shù)模型[17],廣泛應(yīng)用于金融和經(jīng)濟領(lǐng)域。二元Clayton copula定義為

(2)

式中，β∈(-1,+∞){0}。Clayton copula是常見的單變量copula模型,生成函數(shù)為

滿足φ(C(v1,v2))=φ(v1)+φ(v2)性質(zhì)。

三元Clayton copula[19]定義為

(3)

式中,β>0是唯一參數(shù),表示任意2種變量(v1和v2;v1和v3;v2和v3)的相關(guān)性。該模型有很強的對稱性并且結(jié)構(gòu)簡單,適用于三元數(shù)據(jù)建模。

2 基于copula的三元Poisson模型

實際生產(chǎn)過程中,多種質(zhì)量特性在同一時刻的采樣值往往存在相關(guān)信息。某一時刻t采集到M個離散的觀測值,記做Yt=(Yt,1,…,Yt,M)。假設(shè)Yt服從多元Poisson過程,不同時刻采集的數(shù)據(jù)間相互獨立(如Yt,Yt-1),且每個時刻的觀測值存在空間相關(guān)性(如Yt,1,Yt,2)?？紤]Yt的空間相關(guān)性,用copula來構(gòu)建其聯(lián)合概率質(zhì)量函數(shù)。假設(shè)對所有t=1,…,T,向量Yt=(Yt,1,…,Yt,M)的聯(lián)合分布為

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

(4)

式中：β代表copula參數(shù)；Fi代表邊緣累積分布函數(shù),i=1,…,M。

對于二元Poisson過程,在時刻t獲得的2個觀測值可以表達為向量Yt=(Yt,1,Yt,2),其聯(lián)合分布和概率質(zhì)量函數(shù)分別表示為

(5)

同理,三元Poisson過程中變量Yt=(Yt,1,Yt,2,Yt,3)的聯(lián)合分布和概率質(zhì)量函數(shù)分別表示為

(6)

三元Poisson過程的協(xié)方差矩陣表示為

(7)

式中，apq=aqp且apq=cov(Y*,p,Y*,q),p,q=1,2,3。根據(jù)協(xié)方差的定義結(jié)合copula理論可以計算

式中，E(Y*,p)是Y*,p邊緣分布的均值,(Y*,p,Y*,q)的二元聯(lián)合概率質(zhì)量函數(shù)為

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

顯然,對于Poisson分布每個維度變量的取值范圍是非負整數(shù)集。由于各變量的協(xié)方差不恒定,本模型也屬于異方差模型。

實踐中,分布參數(shù)λ及copula參數(shù)β往往是未知的,需要根據(jù)一系列穩(wěn)態(tài)觀測值估計這些參數(shù)。兩階段最大似然估計法,也稱作IFM(inference functions for margins)方法[20],是一種常用的多元copula模型的參數(shù)估計方法,表述如下:

(8)

(9)

兩階段估計法在連續(xù)情形效果顯著,并得到了廣泛應(yīng)用。而在離散情形下,估計參數(shù)在小維度情形很有效,對于大維度數(shù)據(jù)的效果次之[20]。

3 三元CUSUM控制圖設(shè)計

針對多變量監(jiān)控問題,聯(lián)合監(jiān)控多個變量比單獨監(jiān)控一個變量更加容易發(fā)現(xiàn)過程異常[19]。考慮到空間相關(guān)性建模并獲得適用的受控模型及參數(shù)后,需要結(jié)合數(shù)據(jù)特點及控制圖原理設(shè)計出適用的控制圖。本文結(jié)合第2節(jié)的三元相關(guān)Poisson數(shù)據(jù)模型,設(shè)計多元CUSUM控制圖。假設(shè)同一時刻的三元Poisson數(shù)據(jù)存在空間相關(guān),而不同時刻的數(shù)據(jù)間相互獨立。記t時刻的三元Poisson數(shù)據(jù)為Yt=(Yt,1,Yt,2,Yt,3),t=1,…,T,各個維度的參數(shù)記做λ=(λ1,λ2,λ3)。假設(shè)受控狀態(tài)的Poisson參數(shù)為λ0，偏移發(fā)生后變化為λc，變點模型可表示為

根據(jù)似然比檢驗的思想,對變點存在和不存在2種情形的聯(lián)合密度函數(shù)做似然比,構(gòu)造統(tǒng)計量

(10)

式中,f(Ym,…,Yn;λ0)和f(Ym,…,Yn;λc)是(Ym,…,Yn)分別在λ=λ0和λ=λc時的聯(lián)合密度函數(shù)。記

(11)

式中，f(Yi)是概率質(zhì)量函數(shù)。對于本文所討論的二元Poisson過程,f(Yi)=hB(yi,1,yi,2);同理,三元過程中,f(Yi)=hT(yi,1,yi,2,yi,3)?；诖?似然比統(tǒng)計量可以進一步寫做

(12)

等價的遞歸形式可表示為

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

(13)

式中，S0=0。一旦Sn超出給定的控制限h則會報警,意味著存在變點m,即Poisson參數(shù)λ發(fā)生了改變;否則,認為過程參數(shù)λ沒有顯著改變,即過程處于受控狀態(tài)。

4 仿真實驗與分析

4.1 多元相關(guān)Poisson數(shù)據(jù)的生成

基于copula模型生成多元數(shù)據(jù)是建立在條件分布基礎(chǔ)上的。由于copula模型的多樣性,數(shù)據(jù)生成存在一定差異。本節(jié)以三元Poisson過程為例,生成基于Clayton copula的隨機數(shù)據(jù)。記變量v1和(v1,v2)的條件概率密度函數(shù)分別為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數(shù)據(jù)的過程如下所示：

1) 從均勻分布U(0,1)中隨機生成v1;

2) 從均勻分布U(0,1)生成隨機數(shù)作為C2(v2;v1)的實際值,結(jié)合v1求解方程得v2;

3) 從均勻分布U(0,1)生成隨機數(shù)作為C3(v3;v1,v2)的實際值,結(jié)合v1和v2求解方程得v3;

4) 應(yīng)用Poisson分布的反函數(shù)得到一組數(shù)據(jù)(y1,y2,y3),其中,y1=F-1(vi),i=1,2,3。

4.2 參數(shù)估計

4.3 馬爾科夫鏈近似求解平均運行鏈長

目前,主流的用于評價控制圖監(jiān)控效果的指標是平均運行鏈長(average run length,ARL),即運行鏈長(RL)的期望值。而運行鏈長RL是指從開始到第一次失控報警的樣本點數(shù)目。本文采用馬爾科夫鏈(Markov chain)近似計算多元CUSUM控制圖平均運行鏈長。控制限記做h,受控的多元Poisson參數(shù)為λ0=(λ1,λ2,λ3),目標參數(shù)記做λc=θλ0,失控參數(shù)λd=(δ1λ1,δ2λ2,δ3λ3)。首先將受控區(qū)間分為ne等分,每一段的長度Δ=h/ne,每部分區(qū)間可表示為[(k-1)Δ,kΔ),k=1,2,…,ne。一旦統(tǒng)計量落入?yún)^(qū)間k,認為統(tǒng)計量近似取區(qū)間中點值,即(k-0.5)Δ。假定轉(zhuǎn)移矩陣P中的每一個元素pjk代表統(tǒng)計量從區(qū)間j轉(zhuǎn)移到區(qū)間k的概率,記做

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

(17)

式中，j,k=1,2,…,ne。轉(zhuǎn)移概率的計算首先需要列舉出所有可能的三元Poisson數(shù)據(jù)組合,并分別計算其聯(lián)合概率。需要注意的是,計算受控ARL(L0)時,采用hT(Yt;λ0,β);反之,計算失控ARL(L1)時,采用hT(Yt;λd,β)計算。對轉(zhuǎn)移矩陣做一些變換即可得到表示ARL的向量L

L=(I-P)-11

(18)

式中,I是單位矩陣,1是全為1的列向量。由于初始值S0=0,因此向量L的第一個元素即為ARL的近似值。值得注意的是,該方法近似計算ARL的精確程度受到間隔的等分數(shù)ne影響。ne越大,計算的ARL越接近真實值,但是計算時間會相應(yīng)增加。

4.4 三元CUSUM控制圖的監(jiān)控效果

在缺陷監(jiān)測和疾病預(yù)防的背景下,本文關(guān)注多元Poisson數(shù)據(jù)的異常增長,采用單側(cè)向上的控制圖考察監(jiān)控效果。同理,單側(cè)向下或雙側(cè)控制圖可以根據(jù)實際情況選擇。鑒于多元Poisson數(shù)據(jù)的特點,控制圖的控制下限不做要求,控制上限記做h。在階段一,給定適當?shù)氖芸谹RL, 采用馬爾科夫鏈的方法搜索控制限h。

對三元Poisson過程的每個維度分別考慮3種異常參數(shù),則會有33=27種組合。顯然,實現(xiàn)每一種組合的異常參數(shù)集合是不現(xiàn)實的。因此,考慮了9種組合方式。每一種組合方式中,參數(shù)變化幅度可記為δ=(δ1,δ2,δ3),又因為受控參數(shù)λ0=(λ1,λ2,λ3),則異常參數(shù)λd=(δ1λ1,δ2λ2,δ3λ3)。此外,控制圖設(shè)計中的目標參數(shù)λc=(θλ1,θλ2,θλ3)。考慮到設(shè)計過程的便利性,設(shè)定為1.1,1.5和2.0。

表1 L0=200時多元CUSUM控制圖在3種copula參數(shù)和3種目標參數(shù)下的控制限

設(shè)定λ0=(10,10,10)且L0=200,表1和表2分別展示了3種copula參數(shù)和3種目標參數(shù)條件下的控制限和失控平均運行鏈長。其中,序號1～3僅有1個維度的參數(shù)發(fā)生向上階躍,4～6有2個維度發(fā)生變化,7～9表示3個維度均發(fā)生變化的情形。文獻[4,7]采用一種表征馬氏距離的變化尺度(shift size)量,或稱為非中心性參數(shù)(non-centrality para-meter)

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

式中,λ0和λc分別是受控和失控的參數(shù)向量。如公式(7)所示,Σ是協(xié)方差矩陣。非中心性參數(shù)能夠度量過程參數(shù)的整體變化。對于2組異常參數(shù),如果其非中心性參數(shù)相等,則L1也相當。例如,表2中序號1代表的異常參數(shù)的變化幅度為δ=(1.1,1,1),其性能與δ=(1,1.1,1)及δ=(1,1,1.1)情形相當。 同理,序號6代表的異常參數(shù)的變化幅度為δ=(2,1,2),其性能與δ=(2,2,1)及δ=(1,2,2)情形相當。因此,表2所列出的9種組合方式足夠代表全部27種異常情形的ARL性能。從表2中可以看出,當參數(shù)有變化時,L1值均小于L0值,即L0=200,說明本文方法可以檢測多元Poisson過程均值向上偏移。并且,隨著參數(shù)變化的幅度值不斷增大,L1越來越小。對于一種特定的失控狀態(tài),隨著copula參數(shù)的增加,平均運行鏈長普遍增大。也就是說,多元過程的空間相關(guān)性越來越強時,監(jiān)控異常變化的難度也越來越大。對比3種目標幅度值來看,θ=1.1的控制圖對δ1=δ2=δ3=1.1(序號7)的異常情況效果好;θ=1.5的控制圖對δ1=δ2=δ3=1.5(序號8)效果好;θ=2的控制圖對δ1=δ2=δ3=2(序號9)效果好。因此,設(shè)計多元Poisson控制圖時,對較小的變化幅度需要選擇較小的θ, 反之,對于可能存在的較大變化,選擇較大的θ可獲得更理想的效果。

表2 多元CUSUM控制圖在3種copula參數(shù)和3種目標參數(shù)下的失控平均運行鏈長

圖1展示了Poisson數(shù)據(jù)及其在多元CUSUM控制圖的應(yīng)用。設(shè)定異常參數(shù)的變化幅度δ=(1.1,1.1,1.1),則目標參數(shù)的變化幅度選為θ=δ。三元Poisson數(shù)據(jù)由4.1節(jié)方法生成。圖1a)～1c)分別代表β=1,3,5時對應(yīng)的多元Poisson數(shù)據(jù),而對應(yīng)的CUSUM統(tǒng)計量繪制在圖1d)～1f)中。每組Poisson數(shù)據(jù)包括3個維度,每個維度有20個受控數(shù)據(jù)和20個失控數(shù)據(jù),圖中用藍色虛線分隔。紅色水平線代表在L0=200水平下的控制限。從表1中可得出在β分別取1,3,5時對應(yīng)的控制限為2.5, 2.49和2.55。從圖1中可以看出,對于所列出的3種β取值,提出的方法對失控數(shù)據(jù)都有較好的檢出力。

圖1 Poisson過程數(shù)據(jù)及其在CUSUM控制圖的應(yīng)用(L0=200)

圖2對比了本文所提方法與文獻[3]中D圖的性能。為了便于對比,實驗假定δ=δ1=δ2=δ3,其中δ=1代表受控狀態(tài),δ=1.1,1.2,…,2代表失控狀態(tài)。圖2a)～2c)分別展示在β=1,3,5條件下D圖和本文所提多元CUSUM控制圖在θ=1.1,1.5的ARL對比。從圖中可以看出,D圖在β=1條件下監(jiān)控較大偏移量(θ≥1.2)性能優(yōu)于本文所提CUSUM控制圖。當β=3,5時,本文提出的CUSUM控制圖對所有偏移量均取得更小的L1。其中θ=1.1的CUSUM控制圖對較小的偏移量監(jiān)測效果最好,θ=1.5的CUSUM控制圖對較大的偏移量監(jiān)測效果最好。這一點與表2的結(jié)論一致。

圖2 D圖與本文提出CUSUM控制圖的性能對比

5 結(jié) 論

針對具有空間對稱相關(guān)結(jié)構(gòu)的三元Poisson過程監(jiān)控問題,首先采用三元Clayton copula函數(shù)描述其相關(guān)性,并構(gòu)建聯(lián)合分布;在假設(shè)copula參數(shù)不變的前提下,基于對數(shù)似然比檢驗方法結(jié)合聯(lián)合分布信息,推導(dǎo)出多元CUSUM控制圖的監(jiān)控統(tǒng)計量;在給定初始值的前提下,利用馬爾科夫鏈方法近似求解ARL,驗證其有效性,并與D圖做對比。仿真結(jié)果表明,提出的方法能夠有效監(jiān)測過程均值的向上漂移,并且目標偏移量與實際偏移量相等時可以獲得更好的監(jiān)控性能。與D圖相比,當數(shù)據(jù)間的相關(guān)性較強時,提出的方法對于所有的偏移量能取得更好監(jiān)控效果。在未來工作中,將進一步研究時間和空間關(guān)聯(lián)同時存在的多元數(shù)據(jù)建模及監(jiān)控問題。