基于分位數(shù)回歸方法的函數(shù)型數(shù)據(jù)phaseⅠ控制圖

李晴晴，訾雪旻

（天津職業(yè)技術(shù)師范大學(xué)理學(xué)院，天津 300222）

在工業(yè)過程生產(chǎn)中，質(zhì)量特征不僅僅表現(xiàn)為一個(gè)或多個(gè)獨(dú)立變量，某些情形下則通過變量之間的某種函數(shù)關(guān)系來表示，并把此類數(shù)據(jù)稱為函數(shù)型數(shù)據(jù)。近年來，用函數(shù)型數(shù)據(jù)刻畫產(chǎn)品的某些特征進(jìn)而分析和監(jiān)控此類型的數(shù)據(jù)成為統(tǒng)計(jì)過程控制研究的熱點(diǎn)之一，文獻(xiàn)[1-4]介紹了統(tǒng)計(jì)過程控制的一些成果。目前對(duì)于函數(shù)型數(shù)據(jù)，大多采用普通最小二乘法來估計(jì)回歸系數(shù)，但研究者也越來越關(guān)心解釋變量與響應(yīng)變量分布的中位數(shù)、分位數(shù)呈何種關(guān)系。為了更好地刻畫解釋變量在各個(gè)位置上對(duì)響應(yīng)變量的影響，Koenker等[5]提出了分位數(shù)回歸的思想，這種方法能夠更加全面地描述響應(yīng)變量條件分布的全貌，進(jìn)而受到廣大研究者的關(guān)注，文獻(xiàn)[6-9]給出了分位數(shù)回歸在金融、教育、醫(yī)療衛(wèi)生等領(lǐng)域的應(yīng)用，但分位數(shù)回歸應(yīng)用在統(tǒng)計(jì)過程控制中的研究則甚少。本文針對(duì)上述問題，給出分位數(shù)回歸結(jié)合控制圖的研究。假設(shè)所有的樣本函數(shù)型數(shù)據(jù)都來自一個(gè)線性模型，而備擇假設(shè)為在第m1（m1=1，2，…，m-1）個(gè)樣本之后回歸系數(shù)發(fā)生了漂移。對(duì)于上述假設(shè)提出了基于似然比檢驗(yàn)（likelihood ratio test，LRT）的變點(diǎn)法來監(jiān)控回歸系數(shù)的漂移。許多研究者在回歸建模的背景下研究了變點(diǎn)問題，假設(shè)在任何一個(gè)觀察之后都可能存在1 個(gè)變點(diǎn)，利用似然比法檢測(cè)簡(jiǎn)單線性回歸模型中存在的變點(diǎn)，估計(jì)回歸模型中變點(diǎn)位置，以及估計(jì)變點(diǎn)前后的回歸系數(shù)，相關(guān)研究參見文獻(xiàn)[10-13]。正如Woodall[14]所討論的，第Ⅰ階段控制圖表現(xiàn)的好壞通常根據(jù)出現(xiàn)失控點(diǎn)（在控制線以外的點(diǎn)）的概率來衡量的，而這個(gè)概率通常被稱為“誤報(bào)率”。本文在給定誤報(bào)率α=0.04 的情況下，研究基于分位數(shù)回歸的變點(diǎn)法對(duì)函數(shù)型數(shù)據(jù)的監(jiān)控表現(xiàn)，研究的重點(diǎn)是在函數(shù)型數(shù)據(jù)的基礎(chǔ)上構(gòu)建穩(wěn)健的分位數(shù)回歸估計(jì)，利用分段回歸技術(shù)結(jié)合變點(diǎn)法給出似然比檢驗(yàn)統(tǒng)計(jì)量，對(duì)多個(gè)異常點(diǎn)給出檢驗(yàn)和診斷。

1 第Ⅰ階段的線性函數(shù)型數(shù)據(jù)模型假設(shè)

1.1 單樣本模型假設(shè)

對(duì)于單變量的線性函數(shù)型數(shù)據(jù)，假設(shè)解釋變量X和響應(yīng)變量Y 之間的關(guān)系用以下模型表示：在分段的簡(jiǎn)單線性回歸模型中，假設(shè)數(shù)據(jù)集是由單個(gè)樣本（{x1，y1），（x2，y2），…，（xN，yN）}組成的。解釋變量X 與響應(yīng)變量Y 之間的關(guān)系用以下k 段回歸模型表示：

式中：β0j、β1j分別為第j 個(gè)樣本的截距項(xiàng)和斜率項(xiàng)；θj項(xiàng)為段與段之間的變點(diǎn)（通常θ0=0，θk=N）；εi為隨機(jī)誤差項(xiàng)。

分段線性回歸方法可用于檢測(cè)給定樣本中回歸系數(shù)的變化，估計(jì)變點(diǎn)的位置（θj項(xiàng)），并確定適當(dāng)數(shù)量的變點(diǎn)。Hawkins[15]給出了一般分段多元回歸模型的似然公式，并給出了一種動(dòng)態(tài)規(guī)劃算法，用于確定一般分段多元回歸模型的最大似然統(tǒng)計(jì)量。這種動(dòng)態(tài)規(guī)劃算法既適用于同方差模型，也適用于異方差模型。

1.2 多樣本模型假設(shè)

在具有單個(gè)解釋變量X 和響應(yīng)變量Y 的函數(shù)型數(shù)據(jù)數(shù)據(jù)集中，有m 個(gè)樣本，樣本的形式為（{Xi1，Yi1），i=1，2，…，n1}，（{Xi2，Yi2），i=1，2，…，n2}，…，（{Xim，Yim），i=1，2，…，nj}，nj＞2，j=1，2，…，m 即相當(dāng)于第j個(gè)樣本中含有nj個(gè)觀測(cè)點(diǎn)。解釋變量X 與響應(yīng)變量Y之間的關(guān)系為：

式中：εij為隨機(jī)誤差項(xiàng)。

在這種情況下，假設(shè)每個(gè)樣本中沒有系數(shù)變化，重點(diǎn)是檢測(cè)回歸系數(shù)從一個(gè)樣本到另一個(gè)樣本的變化。為了監(jiān)控函數(shù)型數(shù)據(jù)數(shù)據(jù)集中回歸系數(shù)的變化，把所有的樣本合并成一個(gè)樣本容量為的樣本，并對(duì)該樣本在模型（1）的基礎(chǔ)上做分段回歸。正如之前假設(shè)的在每個(gè)樣本間沒有系數(shù)發(fā)生變化，因此模型（1）中的θj項(xiàng)受指示變量i 的限制。針對(duì)上述問題，使用文獻(xiàn)[16]中的LRT 檢驗(yàn)回歸系數(shù)中出現(xiàn)的變化，并且可以遞歸地識(shí)別數(shù)據(jù)集中的多個(gè)變化。

2 估計(jì)樣本系數(shù)的回歸方法

普通的最小二乘回歸方法是線性回歸模型中最基本、最經(jīng)典的方法。它刻畫了響應(yīng)變量的均值，給出了響應(yīng)變量和解釋變量之間的關(guān)系，但有時(shí)人們更關(guān)心不同分位數(shù)下響應(yīng)變量和解釋變量的關(guān)系。相比較普通的最小二乘法，分位數(shù)回歸能夠更加全面地描述響應(yīng)變量條件分布的全貌，也可以分析響應(yīng)變量的中位數(shù)以及其他分位數(shù)如何被解釋變量影響。

前文給出單樣本和多樣本的回歸模型，對(duì)于單樣本模型直接對(duì)其分段，用分位數(shù)回歸的方法估計(jì)其回歸系數(shù)，而對(duì)于多樣本的模型將其合并成一個(gè)樣本容量為N 的樣本對(duì)其進(jìn)行分段，方法同上。正如之前假設(shè)的線性模型，為了便于計(jì)算和表達(dá)，無論是單樣本還是多樣本數(shù)據(jù)集都將其轉(zhuǎn)化為以下形式求其回歸系數(shù)：

設(shè)隨機(jī)變量的分布函數(shù)為F（y）=P（Y≤y），對(duì)于任意的0 ＜τ ＜1 有F-1（τ）=inf{y：F（y）≥τ}，稱τ 為y的分位數(shù)。對(duì)于回歸而言，就是使函數(shù)型數(shù)據(jù)值與真值之間距離最短，隨機(jī)變量Y=（y1，y2，…，yn）的樣本均值回歸使殘差平方和最小，即：

樣本中位數(shù)回歸使殘差絕對(duì)值之和達(dá)到最小，即：

而對(duì)于樣本的分位數(shù)回歸使加權(quán)的殘差和達(dá)到最小，即：

上式可以等價(jià)表示為：

式中：ρτ（u）為檢驗(yàn)函數(shù)，定義為ρτ（u）=u（τ-（Iu ＜0））。

對(duì)于條件分位數(shù)，可以將條件分位數(shù)表示為Qy（τ|x）=xTβ（τ）并且可以通過求解來得到（τ），即：

對(duì)于模型（2）中回歸系數(shù)的解為：

3 基于分位數(shù)回歸的變點(diǎn)控制圖

觀測(cè)數(shù)據(jù)集為m 個(gè)隨機(jī)樣本，每個(gè)樣本由nj對(duì)觀測(cè)序列（xij，yi）j組成，i=1，2，…，n；j=1，2，…，m。解釋變量X 和響應(yīng)變量Y 之間的關(guān)系用模型（2）來表示，并假設(shè)模型中εij項(xiàng)是獨(dú)立同分布于N（0，1）的隨機(jī)變量，進(jìn)一步假設(shè)模型中的設(shè)計(jì)陣X 已知，且每個(gè)樣本中的x 值相同。如果過程是穩(wěn)定的，那么回歸系數(shù)滿足為了監(jiān)測(cè)線性函數(shù)型數(shù)據(jù)集中的回歸系數(shù)變化，將m個(gè)樣本所有觀測(cè)值合并成一個(gè)容量為N 的樣本，然后用對(duì)合并后的樣本進(jìn)行分段回歸。如果樣本中存在多變點(diǎn)，那么當(dāng)監(jiān)測(cè)到一個(gè)變點(diǎn)后，可以用二分法來繼續(xù)尋找其他的變點(diǎn)，直到?jīng)]有變點(diǎn)出現(xiàn)。假設(shè)從第m1（m1=1，2，…，m-1）個(gè)函數(shù)型數(shù)據(jù)之后回歸模型系數(shù)中存在變點(diǎn)，那么有以下假設(shè)。

假設(shè)1

H1：存在一個(gè)（jj=1，2，…，m-1），使得

假設(shè)2

H1：存在一個(gè)（jj=1，2，…，m-1），使得

為了檢驗(yàn)上述假設(shè)，構(gòu)造以下似然比檢驗(yàn)統(tǒng)計(jì)量：

可控狀態(tài)下，式（5）中的lrt m1（τ）統(tǒng)計(jì)量在m1取值不同時(shí)，其期望值也不同，模擬可控狀態(tài)下τ 取0.5和0.9 時(shí)的lrt m1（τ）統(tǒng)計(jì)量的期望值如表1 所示。模擬的可控模型中系數(shù)取值為β0=0，β1=1，即模型為yij=xi+εij（i=1，2，…，n；j=1，2，…，m）。其中n=10，m=20，解釋變量X 的取值為0（0.2）1.8（即0～1.8，間隔為0.2）。假設(shè)εij是獨(dú)立同分布于N（0，1）的隨機(jī)變量。

表1 模擬可控狀態(tài)τ 下取0.5 和0.9 時(shí)的lrt m1（τ）統(tǒng)計(jì)量的期望值

正如表1 所示，m1取值不同時(shí)E（lrt m1（τ））的值也不同，用lrt m1（τ）除以歸一化因子Cm1（τ）矯正lrt m1（τ）使得對(duì)于不同的m1值似然比檢驗(yàn)統(tǒng)計(jì)量的期望值是一致的。如果令Cm1（τ）=E（lrt m1（τ）），則修正后的統(tǒng)計(jì)量為：

然后在第I 類錯(cuò)誤概率α 給定的情況下，通過模擬近似給出修正后的lrt 統(tǒng)計(jì)量的閾值。

4 模擬研究控制圖表現(xiàn)

使用R 軟件進(jìn)行模擬，可控模型中系數(shù)取值為β0=0，β1=1，即模型為yij=xi+εij，其中i=1，2，…，n；j=1，2，…，m；n=10；m=20。假設(shè)εij是獨(dú)立同分布于N（0，1）的隨機(jī)變量。解釋變量為固定樣本，這里假設(shè)解釋變量X 的取值為0（0.2）1.8（即0～1.8，間隔為0.2）。取分位數(shù)τ 為0.5 和0.9，使用式（6）中的似然比檢驗(yàn)統(tǒng)計(jì)量通過10 000 次模擬在第Ⅰ類錯(cuò)誤概率為0.04 的情況下，產(chǎn)生的閾值分別為7.36 和22.02。在模擬中所研究的漂移類型是在樣本q（q ＜m）之后發(fā)生的持續(xù)性漂移，假設(shè)漂移發(fā)生在樣本q=15，18 和19之后，截距和斜率項(xiàng)分別以為單位，其中將加了漂移的似然比檢驗(yàn)統(tǒng)計(jì)量與閾值比較，如果大于閾值，就記為一次成功，運(yùn)行上述程序10 000 次，能準(zhǔn)確識(shí)別變化的概率稱為整體失控報(bào)警概率。

截距從β0漂移到的失控報(bào)警概率如圖1 所示。圖1（a）和（b）分別給出了τ 取0.5 和0.9時(shí)截距在樣本q=15、18 和19 之后（q ＜m）從β0漂移到情況下的整體失控報(bào)警概率，其中λ取0.5（0.5）5（即0.5～5 間隔為0.5）。從圖1 可以看出，本文提出的變點(diǎn)法在q 取某一值的情況下，監(jiān)控大的漂移比監(jiān)控小的漂移好，在漂移取固定值時(shí)q 取值越小變點(diǎn)法的監(jiān)控表現(xiàn)就越好。但當(dāng)q 值取19 時(shí)，從圖1 中可以看出變點(diǎn)法的監(jiān)控效果明顯弱于q 取其他值的監(jiān)控效果。

圖1 截距從β0 漂移到β0+λσ/的失控報(bào)警概率

斜率項(xiàng)從β1漂移到的失控報(bào)警概率如圖2 所示。圖2（a）和（b）分別給出了取0.5 和0.9時(shí)斜率在樣本q=15、18 和19 之后（q ＜m）從β1漂移到情況下的整體失控報(bào)警概率，其中δ取從圖2 可以看出，對(duì)于監(jiān)控這種類型的漂移結(jié)果與監(jiān)控截距項(xiàng)結(jié)果一致，都是在q 取某一值的情況下，監(jiān)控大的漂移比監(jiān)控小的漂移好，在漂移δ 取固定值時(shí)，q 取值越小變點(diǎn)法的監(jiān)控表現(xiàn)就越好，但是當(dāng)δ 取值足夠大時(shí)，即截距項(xiàng)發(fā)生大漂移時(shí)，無論q 取何值，變點(diǎn)法的監(jiān)控結(jié)果都令人滿意。

圖2 斜率項(xiàng)從β1 漂移到β1+δσ/的失控報(bào)警概率

截距項(xiàng)從β0漂移到以及斜率項(xiàng)從β1漂移到的失控報(bào)警概率如圖3 所示。圖中給出了τ 取0.5 時(shí)在樣本q=2 和10，截距項(xiàng)從β0漂移到情況下的整體失控報(bào)警概率和斜率項(xiàng)從β1漂移到的失控報(bào)警概率。從圖3 中可以看出，變點(diǎn)控制圖監(jiān)控截距漂移過程并不理想，在q 值取2 或者更大值時(shí)，依然無法準(zhǔn)確監(jiān)測(cè)到變點(diǎn)，而對(duì)于斜率項(xiàng)的漂移來說，變點(diǎn)法可以監(jiān)測(cè)出變點(diǎn)，但效果不如監(jiān)控持續(xù)性漂移好。

由圖1 和圖2 可知，本文提出的基于分位數(shù)回歸的變點(diǎn)法控制圖可以有效監(jiān)控系統(tǒng)變化，并且可以說明使用分位數(shù)回歸法估計(jì)回歸系數(shù)的穩(wěn)健性。

圖3 截距項(xiàng)從β0 漂移到β0+λσ/以及斜率項(xiàng)從β1 漂移到β1+δσ/的失控報(bào)警概率

5 結(jié) 論

本文研究了函數(shù)型數(shù)據(jù)的變點(diǎn)問題，并提出了基于分位數(shù)回歸的變點(diǎn)法，將這種變點(diǎn)法用來監(jiān)控phaseⅠ中系統(tǒng)變化。通過模擬計(jì)算，將監(jiān)控持續(xù)性漂移和非持續(xù)性漂移表現(xiàn)進(jìn)行比較，以及在不同樣本q之后加漂移進(jìn)行比較。結(jié)果表明：無論是持續(xù)性漂移還是非持續(xù)性漂移，q 值越小變點(diǎn)法的表現(xiàn)越好，而本文提出的變點(diǎn)法在監(jiān)控持續(xù)性漂移表現(xiàn)要比監(jiān)控非持續(xù)性漂移表現(xiàn)效果好；當(dāng)回歸系數(shù)發(fā)生較大漂移時(shí)，無論q 取何值、漂移類型是持續(xù)性還是非持續(xù)性，變點(diǎn)法的監(jiān)控結(jié)果都是令人滿意的。