多元正態(tài)向量變點(diǎn)在線監(jiān)測

劉維奇，劉晉芳，秦瑞兵

（1.山西大學(xué) 管理與決策研究所，山西 太原 030006；2.山西大學(xué) 數(shù)學(xué)科學(xué)學(xué)院，山西 太原 030006）

0 引言

在質(zhì)量工程控制中，需要對生產(chǎn)線上的產(chǎn)品進(jìn)行抽檢以檢驗(yàn)產(chǎn)品質(zhì)量是否合格，特別是要檢驗(yàn)產(chǎn)品的質(zhì)量是否發(fā)生顯著波動或超出其質(zhì)量控制范圍，從而當(dāng)產(chǎn)品的質(zhì)量發(fā)生突變時能及時報警，避免出現(xiàn)更多的次品。出現(xiàn)質(zhì)變的時刻稱為變點(diǎn)。在實(shí)際應(yīng)用中，數(shù)據(jù)的生成機(jī)制常由于內(nèi)在或外在因素的變化而發(fā)生變化。因此很多實(shí)際數(shù)據(jù)存在變點(diǎn)，而對變點(diǎn)的統(tǒng)計推斷可以使數(shù)據(jù)的統(tǒng)計分析更為準(zhǔn)確，因此變點(diǎn)從質(zhì)量控制領(lǐng)域逐漸應(yīng)用到經(jīng)濟(jì)、金融、醫(yī)學(xué)、氣候分析等諸多領(lǐng)域。變點(diǎn)問題最初由Page［1］提出，由于其理論和應(yīng)用的重要性而受到很多統(tǒng)計學(xué)者的關(guān)注。Bai［2］基于最小二乘法給出了線性過程均值變點(diǎn)的估計。Bai［3］用極大似然法考慮回歸模型系數(shù)的多變點(diǎn)問題。Perron 基于最小二乘法考慮了對參數(shù)有約束的線性回歸模型的多變點(diǎn)問題。近年來，多元時間序列的變點(diǎn)問題也逐漸受到統(tǒng)計學(xué)者的關(guān)注。Horváth［5］提出了監(jiān)測多元相依穩(wěn)定過程均值變點(diǎn)的幾個統(tǒng)計量。Boutahar［6］提出了一個檢驗(yàn)多元時間序列的非參數(shù)CUSUM方法。關(guān)于變點(diǎn)的其他文獻(xiàn)，參考 Cs?rg?和 Horváth［7］，Perron［8］。

變點(diǎn)問題可以分為兩類：一類是后驗(yàn)檢驗(yàn)，即在給定的一組固定樣本中檢測是否有變點(diǎn)發(fā)生；一類是在線監(jiān)測，即在已有模型的基礎(chǔ)上對新觀測的數(shù)據(jù)監(jiān)測，直到出現(xiàn)變點(diǎn)才停止。在實(shí)際應(yīng)用中對金融數(shù)據(jù)的監(jiān)測更多采用實(shí)時監(jiān)測。此即Chu［9］提出的在線監(jiān)測問題。Carsoule和Franses［10］基于極大似然法對歐洲名義匯率波動的方差變點(diǎn)進(jìn)行監(jiān)測。Andreou等［11］和Zeileisd等［12］用最小二乘法考慮了動態(tài)經(jīng)濟(jì)模型中的變點(diǎn)監(jiān)測問題。Horváth等［13］提出了分別基于殘差累積和與遞歸殘差累積和的統(tǒng)計量監(jiān)測線性回歸模型中變點(diǎn)的兩種監(jiān)測方案。其模擬結(jié)果表明，該方法能夠較快地監(jiān)測到較早發(fā)生的變點(diǎn)，而對于離開始位置較遠(yuǎn)的變點(diǎn)其監(jiān)測的平均運(yùn)行長度較長。陳占壽［14］通過引入窗寬參數(shù)改變監(jiān)測的開始時刻對文獻(xiàn)［13］中方法進(jìn)行改進(jìn)，這樣當(dāng)變點(diǎn)出現(xiàn)的時刻離最初監(jiān)測的位置較遠(yuǎn)時，可以通過滑動監(jiān)測的開始時刻使得變點(diǎn)看起來較早發(fā)生。陳占壽［15］用滑動核加權(quán)比率統(tǒng)計量監(jiān)測帶有趨勢項時間序列的持久性變點(diǎn)。

針對多元時間序列均值變點(diǎn)問題，本文采用陳占壽［14］引進(jìn)窗寬參數(shù)的方法，構(gòu)造帶有窗寬參數(shù)的殘差累積和統(tǒng)計量監(jiān)測多元正態(tài)向量均值變點(diǎn)。模擬表明引入窗寬參數(shù)可以改變監(jiān)測開始的時刻，從而將較晚出現(xiàn)的變點(diǎn)移動到較早的時刻以增大檢驗(yàn)的勢，同時能減小平均運(yùn)行長度，達(dá)到快速監(jiān)測變點(diǎn)的目的。

1 模型

我們考慮下列模型：其中Xi是一個d×1維正態(tài)隨機(jī)向量，μi是一個d×1維向量，ei是d×1維殘差序列，滿足下面假設(shè)條件：

假設(shè)1：誤差序列滿足：

（1）｛ei，i＝1，2，…，n｝獨(dú)立同分布，且有E ei，j＝0，（1≤i＜∞，1≤j≤d）和E‖ei‖ν＜∞，ν＞2.

（2）存在正整數(shù)m，使得對任意的k和l，若1－k＞m，則有σ｛ei，1≤i＜k｝和σ｛ei，l≤i＜∞｝相互獨(dú)立。

假設(shè)2：｛μi，i＝1，2，…，n｝和｛ei，i＝1，2，…，n｝相互獨(dú)立。

假設(shè)3：存在一d維矩陣D使得

假設(shè)4：向量Xi的各分量彼此獨(dú)立，記

假設(shè)5：前m個歷史樣本是無污染的，即

以上假設(shè)1－3是證明監(jiān)測統(tǒng)計量的漸近性質(zhì)的必要條件，與Horváth［13］中的假設(shè)相同。

現(xiàn)在基于前m個樣本觀測值，從第m＋1個觀測值開始監(jiān)測，即檢驗(yàn)原假設(shè)（H0）

和備擇假設(shè)（H1）

其中向量μ0，μ1和k＊≥1均未知。

2 監(jiān)測方法和主要結(jié)果

我們構(gòu)造的殘差CUSUM統(tǒng)計量如下：

3 模擬計算和實(shí)例分析

3.1 模擬計算

這部分通過蒙特卡羅模擬方法檢驗(yàn)本文方法有限樣本下的性質(zhì)。因?yàn)樵谀M中矩陣D是未知參數(shù)矩陣，需要估計。由Horváth［5］中結(jié)果知，可用基于前m個樣本觀測值的Dm估計，其中

且滿足：‖Dm－D‖＝oP（（log logm）－1／2）.可以證明前面的定理在Dm下依然成立。

表1 臨界值表Table 1 Critical values

表2 監(jiān)測樣本量q＝2 m時的經(jīng)驗(yàn)水平Table 2 Empirical size when q＝2m

表1列出了當(dāng)d＝3時，基于的50 000次模擬得到的臨界值。表2是樣本量m＝200，監(jiān)測長度q＝2m時的經(jīng)驗(yàn)水平。從表中得知，隨著h的增大，經(jīng)驗(yàn)水平越來越大，也就是說h越大監(jiān)測統(tǒng)計量越靈敏。另外，當(dāng)γ＝0.45，h≥0.2和γ＝0.25，h≥0.4時存在水平失真問題。而在其他情況下，隨著γ和h的增大，經(jīng)驗(yàn)水平逐漸靠近檢驗(yàn)水平。

表3和表4是考慮均值向量從（0，0，0）′變化到（0.6，0.6，0.6）′時模擬得到的勢，其中樣本容量m＝50，100，200監(jiān)測長度q＝m。從表中可得：隨著樣本量的增加，檢驗(yàn)的勢逐漸增大；在樣本量和變點(diǎn)位置都相同的情況下，γ的增大會提高檢驗(yàn)的勢，同時減小平均運(yùn)行長度；當(dāng)變點(diǎn)位置靠前時，在樣本量較大的時候，改變h的取值對檢驗(yàn)的勢和平均運(yùn)行長度影響不大，甚至較大的h值會減小檢驗(yàn)的勢，同時平均運(yùn)行長度會增加；而當(dāng)變點(diǎn)的位置靠后時，引入窗寬參數(shù)會明顯提升檢驗(yàn)的勢，且此時檢驗(yàn)的平均運(yùn)行長度明顯減小，因此本文方法可以適用于變點(diǎn)位置靠后的情況。

表3 檢驗(yàn)的勢（顯著水平為5%）Table 3 Empirical power of the procedure

表4 檢驗(yàn)的平均運(yùn)行長度（顯著水平為5%）Table 4 ARLs of the procedures

此外γ＝0.45，h＝0.2適用于監(jiān)測不久就出現(xiàn)變點(diǎn)的情況，γ＝0.25，h＝0.5適用于變點(diǎn)較晚出現(xiàn)的情況。而在實(shí)際應(yīng)用中，當(dāng)我們不能確定變點(diǎn)的大致位置時，可以選擇γ＝0.25，h＝0.2，以便達(dá)到兩者的折中且盡可能快速地監(jiān)測到變點(diǎn)。

3.2 實(shí)例分析

通過一個白酒的生產(chǎn)過程實(shí)例來說明本文方法的有效性。數(shù)據(jù)包括從2004年5月到2007年2月的4 898個觀測，該數(shù)據(jù)集可以從UCI機(jī)器學(xué)習(xí)數(shù)據(jù)庫（http ∥archive.ics.uci.edu ml datasets Wine＋Quality）中的“Wine Quality Data Set”中得到。此數(shù)據(jù)集中包括12個變量，分別是：非揮發(fā)性酸、揮發(fā)性酸度、檸檬酸、殘留糖、氯化物、游離二氧化硫、總二氧化硫、密度、p H、硫酸鹽和乙醇（分別記作y1，y2，y3，y4，y5，y6，y7，y8，y9，y10，y11）。還有一個由感官分析得到的質(zhì)量分類變量，此分類變量的取值為0－10，取值為6時表示口感最好，取值越偏離6表示口感越差，取值為5－7都表示質(zhì)量在可控范圍之內(nèi)。關(guān)于此數(shù)據(jù)集更多詳細(xì)的信息可以參考Cortez［16］。

我們的目標(biāo)是監(jiān)測VinhoVerde白葡萄酒的生產(chǎn)過程以確保所生產(chǎn)酒的質(zhì)量（口感）合格。為了方便，在這里我們只對前1 000個數(shù)據(jù)做分析，選取三個分量，分別是：非揮發(fā)性酸（y1）、總二氧化硫（y7）、密度（y9）。圖1是這三個變量的原始數(shù)據(jù)圖，圖2的Q-Q圖表明這三個變量分別都服從正態(tài)分布，通過計算相關(guān)系數(shù)矩陣得知此三個變量是不相關(guān)的。圖3的三元Q-Q圖表明該數(shù)據(jù)服從多元正態(tài)分布。因此，我們所選用的數(shù)據(jù)集符合前文中提出的假設(shè)。

Fig.1 Raw data of y1，y7 and y 9圖1 變量y 1、y 7和y 9的原始數(shù)據(jù)散點(diǎn)圖

Fig.2 Q-Q plot of y 1，y 7 and y 9圖2 變量y 1、y 7和y 9的Q-Q圖

Fig.3 Q-Q plot of y 1，y 7 and y 9圖3 變量的三元Q-Q圖

下面開始對該數(shù)據(jù)集做變點(diǎn)監(jiān)測。如圖1所示，前150個數(shù)據(jù)各分量均值中不存在明顯的變化，因此選取前150個數(shù)據(jù)作為歷史樣本。我們把監(jiān)測到的變點(diǎn)位置標(biāo)注到質(zhì)量分類變量的散點(diǎn)圖中，如圖4所示，令γ＝0.25，h＝0，在5%的顯著性水平下從第151個樣本（圖4中紅虛線的位置）開始監(jiān)測。在k＝187的位置（圖4中紅實(shí)線的位置）監(jiān)測到變點(diǎn)。若選取γ＝0.25，h＝0.25，則在k＝154的位置（圖4中綠實(shí)線的位置）監(jiān)測到變點(diǎn)，該結(jié)果表明采用引入窗寬參數(shù)的方法確實(shí)能夠更早的監(jiān)測到變點(diǎn)。

4 結(jié)論與展望

Fig.4 Q-Q plot of y 1、y 7 and y 9圖4 監(jiān)測結(jié)果圖

對于存在變點(diǎn)的多元時間序列，如何較早的監(jiān)測到變點(diǎn)，對控制風(fēng)險做出有效決策具有現(xiàn)實(shí)意義。本文通過采用引入窗寬參數(shù)法，使得多元時間序列均值變點(diǎn)靠后時能盡快監(jiān)測到變點(diǎn)。模擬表明窗寬參數(shù)選取合理可以使檢驗(yàn)的誤報率減小，同時使得平均運(yùn)行長度減小。另外，本文考慮的是分量之間獨(dú)立的變點(diǎn)監(jiān)測，對于分量相依的情形還有待后續(xù)研究。

［1］ Page E S.A Test for a Change in a Parameter Occurring at an Unknown Point［J］.Biometrika，1955，42（3）：523-527.

［2］ Bai J.Least Squares Estimation of a Shift in Linear Processes［J］.Journal of Time Series Analysis，1994，15（5）：453-472.

［3］ Bai J.Likelihood Ratio Tests for Multiple Structural Changes［J］.Journal of Econometrics，1999，91（2）：299-323.

［4］ Perron P，Qu Z.Estimating Restricted Structural Change Models［J］.Journal of Econometrics，2006，134（2）：373-399.

［5］ Horváth L，Kokoszka P，Steinebach J.Testing for Changes in Multivariate Dependent Observations with an Application to Temperature Changes［J］.Journal of Multivariate Analysis，1999，68（1）：96-119.

［6］ Boutahar M.Testing for Change in Mean of Independent Multivariate Observations with Time Varying Covariance［J］.Journal of Probability and Statistics，2012，2012：1-17.

［7］ Cs?rg?M，Horváth L.Limit Theorems in Change-Point Analysis［M］.New York：John Wiley ＆ Sons，1997.

［8］ Perron P.Dealing with Structural Breaks［M］.Palgrave Handbook of Econometrics Vol 1，New York：Palgrave Macmillan，2006：278-352.

［9］ Chu C S J，Stinchcombe M，White H.Monitoring Structural Change［J］.Econometrica，1996，64（5）：1045-1063.

［10］ Carsoule F，F(xiàn)ranses P H.Monitoring Structural Change in Variance，with an Application to European Nomianl Exchange Rate Volatility［J］.Econometric，1999，6：1-33.

［11］ Andreou E，Ghysels E.Monitoring Disruptions in Finacial Markets［J］.Journal of Econometrics，2006，135（1）：77-124.

［12］ Zeileis A，Leisch F，Kleiber C，et al.Monitoring Structural Change in Dynamic Econometric Models［J］.Journal of Applied Econometrics，2005，20（1）：99-121.

［13］ Horváth L，Hu?kováM，Kokoszka P，et al.Monitoring Changes in Linear Models［J］.Journal of Statistical Planning and Inference，2004，126（1）：225-251.

［14］ 陳占壽，田錚，丁明濤.線性回歸模型參數(shù)變點(diǎn)的在線監(jiān)測［J］.系統(tǒng)工程理論與實(shí)踐，2010，30（6）：1047-1054.

［15］ 陳占壽，田錚.含趨勢項時間序列持久性變點(diǎn)監(jiān)測［J］.系統(tǒng)工程理論與實(shí)踐，2014，34（4）：936-943.

［16］ Cortez P，Cerdeira A，Almeida F，et al.Modeling Wine Preferences by Data Mining from Physicochemical Properties［J］.Decision Support Systems，2009，47（4）：547-533.