李拂曉, 肖燕婷, 陳占?jí)?/p>

(1.西安理工大學(xué)理學(xué)院, 陜西 西安 710054;2.青海師范大學(xué)數(shù)學(xué)與統(tǒng)計(jì)學(xué)院,青海 西寧 810008)

1.引言

分類變量在現(xiàn)實(shí)生活中十分常見, 如每個(gè)學(xué)生的期末考試成績(jī)按分?jǐn)?shù)段可分為“優(yōu)秀”,“良好”, “及格”和“不及格”四類; 每個(gè)人的血壓值可分為低血壓, 正常和高血壓三類.按照時(shí)間順序記錄分類變量值, 即構(gòu)成分類時(shí)間序列.關(guān)于分類時(shí)間序列的建模, 研究者提出了多種模型, 如Markov鏈模型, 廣義線性模型, 整數(shù)自回歸模型和自回歸滑動(dòng)平均模型等.廣義線性模型可以把響應(yīng)變量的條件期望表示為過(guò)去觀測(cè)值和協(xié)變量的函數(shù), 且相比Markov鏈模型, 既不需要假設(shè)馬氏性, 也不需要假設(shè)平穩(wěn)性, 是分類時(shí)間序列建模比較常用的模型[1].累積Logistic回歸模型(Cumulative logistic regression model)是一類重要的廣義線性模型, 在生物醫(yī)學(xué), 社會(huì)科學(xué)和遺傳學(xué)等領(lǐng)域有著廣泛的應(yīng)用.

在數(shù)據(jù)分析中，數(shù)據(jù)所服從的模型由于受到某些因素的影響, 可能在某個(gè)時(shí)刻或位置發(fā)生改變, 即產(chǎn)生變點(diǎn).過(guò)去的幾十年中，有許多學(xué)者研究了模型的結(jié)構(gòu)變點(diǎn)檢驗(yàn)和估計(jì)問(wèn)題.[2-9]變點(diǎn)檢驗(yàn)一般分為回顧性檢驗(yàn)(Retrospective test)和序貫檢驗(yàn)(Sequential test), 前者是分析已觀測(cè)的歷史數(shù)據(jù), 檢驗(yàn)數(shù)據(jù)中是否存在變點(diǎn); 后者是連續(xù)檢驗(yàn)數(shù)據(jù)中是否存在變點(diǎn), 即新觀測(cè)的數(shù)據(jù)是否仍符合原有的模型,也稱為在線監(jiān)測(cè)[10].針對(duì)回顧性變點(diǎn)檢驗(yàn),Antoch[11]基于整體最大值型(Overall maximum type)統(tǒng)計(jì)量研究了廣義線性模型的結(jié)構(gòu)變點(diǎn)檢驗(yàn)問(wèn)題.Hudecova[12]利用得分型統(tǒng)計(jì)量研究了二分類自回歸模型(Binary autoregressive model)的結(jié)構(gòu)變點(diǎn)檢驗(yàn)問(wèn)題.Fokianos[13]等基于偏似然得分過(guò)程構(gòu)造統(tǒng)計(jì)量, 檢驗(yàn)二分類Logistic回歸模型的系數(shù)是否存在變點(diǎn), Gombay[14]等將其推廣到多分類Logistic回歸模型.WANG等[15]在高維同質(zhì)性檢驗(yàn)的基:上提出了一種新的方法, 研究了具有大量類別的多分類數(shù)據(jù)的變點(diǎn)估計(jì)問(wèn)題.關(guān)于變點(diǎn)的序貫檢驗(yàn), XIA等[16]分別基于殘差的累計(jì)和(Cumulative sum, CUSUM)和滑動(dòng)和(Moving sum, MOCUSM)構(gòu)造統(tǒng)計(jì)量, 在線監(jiān)測(cè)廣義線性模型的結(jié)構(gòu)是否存在變點(diǎn).Hohle[17]利用CUSUM方法在線監(jiān)測(cè)分類時(shí)間序列所服從的Logistic回歸模型是否存在結(jié)構(gòu)變化.LI等[18]基于偏最大似然得分過(guò)程構(gòu)造統(tǒng)計(jì)量, 研究了多分類Logistic回歸模型的變點(diǎn)序貫檢驗(yàn)問(wèn)題.

本文基于偏似然得分過(guò)程構(gòu)造序貫檢驗(yàn)統(tǒng)計(jì)量, 在線監(jiān)測(cè)累積Logistic回歸模型的結(jié)構(gòu)是否存在變點(diǎn).原假設(shè)下推導(dǎo)出統(tǒng)計(jì)量的漸近分布, 備擇假設(shè)下證明其一致性.與文[17]提出的CUSUM變點(diǎn)監(jiān)測(cè)方法相比, 本文所提出的方法只需要在原假設(shè)下估計(jì)出模型參數(shù), 備擇假設(shè)下不需要重新估計(jì).模擬試驗(yàn)表明原假設(shè)下統(tǒng)計(jì)量的檢驗(yàn)水平大部分都接近于顯著性水平0.05, 在備擇假設(shè)下雖然監(jiān)測(cè)到變點(diǎn)的平均運(yùn)行長(zhǎng)度較長(zhǎng), 但監(jiān)測(cè)到變點(diǎn)的誤報(bào)率較低, 且檢驗(yàn)勢(shì)也較高.最后通過(guò)一組實(shí)際數(shù)據(jù)說(shuō)明本文方法的有效性.

2.模型及變點(diǎn)檢驗(yàn)問(wèn)題

變點(diǎn)的序貫檢驗(yàn)統(tǒng)計(jì)量一般是由監(jiān)測(cè)統(tǒng)計(jì)量Γ(m,k)和邊界函數(shù)g(m,k)構(gòu)成.定義停時(shí)

Γ(m,k)和g(m,k)在原假設(shè)下滿足

其中α ∈(0,1).

3.變點(diǎn)的序貫檢驗(yàn)

首先利用偏最大似然方法對(duì)θ進(jìn)行估計(jì)[20], 偏似然函數(shù)為

則對(duì)數(shù)偏似然函數(shù)為

記

其中

假設(shè)1和假設(shè)2表明偏對(duì)數(shù)似然函數(shù)的二階導(dǎo)數(shù)是關(guān)于θ的連續(xù)函數(shù)[21].假設(shè)2保證了Ut(θ)是非奇異的.假設(shè)3表明

知了在樹上放肆地叫著，夕陽(yáng)的余暉照過(guò)來(lái)，我看到西天堆起大山一樣的烏云，一重連一重。一注一注的陽(yáng)光從烏云的后面射向天空，給綿綿云山綴上了一帶炫目的金邊。我強(qiáng)撐著爬起來(lái)，漫無(wú)目的朝前走著。身后不知么事時(shí)候，冒出幾個(gè)細(xì)伢兒跟著，他們齊聲唱道——

滿足正定性的概率是1.[21]若假設(shè)4成立, 則當(dāng)m →∞時(shí)

其中G(θ)為正定矩陣.[21]

令Sm(θ)=0, 求解可得θ的偏最大似然估計(jì)量基于構(gòu)造變點(diǎn)序貫檢驗(yàn)統(tǒng)計(jì)量

g(s)是邊界函數(shù).

下面兩個(gè)定理給出統(tǒng)計(jì)量Stat在原假設(shè)下的極限分布以及備擇假設(shè)下的一致性.

定理3.1如果假設(shè)1-4和原假設(shè)H0成立, 則

(i) 若g(s)=c, 則

由Nκm」, 類似于文[13]中定理1和定理3的證明可得

則定理的結(jié)論成立.(ii)的證明與(i)類似.

根據(jù)假設(shè)1-4,

因δ0, 當(dāng)m →∞時(shí)

則(i)的結(jié)論成立, 當(dāng)m+1＜k ＜k*時(shí), 證明是類似的.(ii)的結(jié)論可由(i)直接推出.

4.模擬試驗(yàn)

本節(jié)利用Monte Carlo模擬研究變點(diǎn)序貫檢驗(yàn)統(tǒng)計(jì)量Stat的有限樣本性質(zhì).在推導(dǎo)統(tǒng)計(jì)量的漸近分布時(shí), 采用了兩種不同的邊界函數(shù)g(s)=c和g(s)=cg1(s)=cs0.8, 相應(yīng)的臨界值可分別通過(guò)下式獲得：

其中α是顯著性水平.數(shù)據(jù)生成過(guò)程如下：

其中{Yt}是類別數(shù)為3的累積logistic回歸模型,α1=-0.5,α2=0.2,(φ1,φ2,φ3)T=(2,-0.5,1)T.假設(shè)歷史樣本量m= 100,200,500, 監(jiān)測(cè)樣本量N1=N -m= (κ-1)m,κ= 3,5,7.顯著性水平取α=0.05, 模擬試驗(yàn)重復(fù)的次數(shù)為2500 次.

令θ= (α1,α2,φ1,φ2,φ3)T, 在現(xiàn)實(shí)應(yīng)用中可能只關(guān)注部分參數(shù)是否發(fā)生變化.本文假設(shè)只監(jiān)測(cè)參數(shù)α1和φ1是否發(fā)生變化, 其余參數(shù)作為冗余參數(shù).表1給出了原假設(shè)下統(tǒng)計(jì)量Stat在兩種邊界函數(shù)下的檢驗(yàn)水平.結(jié)果表明兩種情形下的檢驗(yàn)水平大部分都接近于顯著性水平0.05, 少部分存在一些扭曲.

表1 兩種邊界函數(shù)下統(tǒng)計(jì)量的檢驗(yàn)水平

考慮下述形式的備擇假設(shè):

HA:α1在m+k*處從-0.5變?yōu)?2,φ1在m+k*處從2變?yōu)?,

其中k*= 0.01N,0.05N,0.1N,0.2N,0.3N.表2-5描述了統(tǒng)計(jì)量Stat的檢驗(yàn)勢(shì)和監(jiān)測(cè)到變點(diǎn)運(yùn)行長(zhǎng)度(Run length)的四分之一分位數(shù)(Q1), 中位數(shù)(Q2), 均值(Average run length, ARL), 四分之三分位數(shù)(Q3), 最大值(Max) 以及誤報(bào)率(Pτ).運(yùn)行長(zhǎng)度是變點(diǎn)發(fā)生位置與監(jiān)測(cè)到變點(diǎn)的位置之間的樣本量, 其均值稱為平均運(yùn)行長(zhǎng)度(ARL).檢驗(yàn)勢(shì)是監(jiān)測(cè)到變點(diǎn)發(fā)生在m+k*之后,N之前的次數(shù)所占的比率.ARL≥ˉm+k*表示監(jiān)測(cè)到變點(diǎn)位置發(fā)生在m+k*之后時(shí), 其運(yùn)行長(zhǎng)度的平均值.誤報(bào)率是指監(jiān)測(cè)到變點(diǎn)發(fā)生在m+k*之前的次數(shù)所占的比率.做為對(duì)比, 考慮了文[17]中CUSUM變點(diǎn)檢測(cè)方法, 其臨界值可通過(guò)模擬平均運(yùn)行長(zhǎng)度的方式獲取.假設(shè)受控狀態(tài)下平均運(yùn)行長(zhǎng)度為370, 則相應(yīng)的臨界值為2.5.CUSUM統(tǒng)計(jì)量的檢驗(yàn)勢(shì), 運(yùn)行長(zhǎng)度的數(shù)字特征和誤報(bào)率見表6.

表2 在HA下, 當(dāng)m=100,N =400, g(s)=c時(shí), Stat的檢驗(yàn)勢(shì), 運(yùn)行長(zhǎng)度的數(shù)字特征和誤報(bào)率

表3 在HA下, 當(dāng)m=100,N =400, g(s)=cs0.8時(shí), Stat的檢驗(yàn)勢(shì), 運(yùn)行長(zhǎng)度的數(shù)字特征和誤報(bào)率

表4 在HA下, 當(dāng)m=200,N =800, g(s)=c時(shí), Stat的檢驗(yàn)勢(shì), 運(yùn)行長(zhǎng)度的數(shù)字特征和誤報(bào)率

表5 在HA下, 當(dāng)m=200,N =800, g(s)=cs0.8時(shí), Stat的檢驗(yàn)勢(shì), 運(yùn)行長(zhǎng)度的數(shù)字特征和誤報(bào)率

從表2-5的結(jié)果可以看出, Stat的檢驗(yàn)勢(shì)隨著樣本量的增加而增大, 但當(dāng)變點(diǎn)發(fā)生的位置離監(jiān)測(cè)起始時(shí)刻越遠(yuǎn), 其檢驗(yàn)勢(shì)越低.邊界函數(shù)不同, 其對(duì)應(yīng)的檢驗(yàn)勢(shì)也有差異.當(dāng)m=100,N=400時(shí), 邊界函數(shù)取g(s)=c的檢驗(yàn)勢(shì)高于g(s)=cs0.8時(shí)的檢驗(yàn)勢(shì).但當(dāng)樣本量增大,m=200,N=800時(shí), 兩種邊界函數(shù)下的檢驗(yàn)勢(shì)都趨于1.

運(yùn)行長(zhǎng)度的數(shù)字特征Q1,Q2, ARL,Q3, Max隨著樣本量的增大而增大.變點(diǎn)發(fā)生的位置也會(huì)對(duì)運(yùn)行長(zhǎng)度產(chǎn)生影響.當(dāng)邊界函數(shù)取g(s) =c時(shí), 運(yùn)行長(zhǎng)度受變點(diǎn)位置的影響較小.例如當(dāng)m=100,N=400, 不同變點(diǎn)位置下的平均運(yùn)行長(zhǎng)度是230左右;m=200,N=800, 不同變點(diǎn)位置下的平均運(yùn)行長(zhǎng)度是293左右.但當(dāng)g(s) =cs0.8時(shí), 不同變點(diǎn)位置對(duì)運(yùn)行長(zhǎng)度的影響較大.例如m=100,N=400,k*=0.01N,0.3N時(shí), 平均運(yùn)行長(zhǎng)度分別是156.66和259.13, 差異較大.相比g(s) =c, 邊界函數(shù)取g(s) =cs0.8時(shí), 若變點(diǎn)的位置離監(jiān)測(cè)起始時(shí)刻較近時(shí), ARL較短, 但當(dāng)變點(diǎn)位置離監(jiān)測(cè)起始時(shí)刻較遠(yuǎn)時(shí), ARL較長(zhǎng).由于誤報(bào)率在大多數(shù)情形下都為0, 變點(diǎn)后的平均運(yùn)行長(zhǎng)度(ARL≥ˉm+k*) 與平均運(yùn)行長(zhǎng)度(ARL)幾乎都是相等的.

表6 在HA下, 當(dāng)m=100,N =400, CUSUM變點(diǎn)檢測(cè)方法的檢驗(yàn)勢(shì), 運(yùn)行長(zhǎng)度的數(shù)字特征和誤報(bào)率

由表6可以看出, 當(dāng)利用CUSUM方法監(jiān)測(cè)變點(diǎn)時(shí), 其檢驗(yàn)勢(shì), 運(yùn)行長(zhǎng)度和誤報(bào)率受變點(diǎn)位置的影響較大.例如m= 100,N= 400, 當(dāng)變點(diǎn)發(fā)生在監(jiān)測(cè)起始時(shí)刻時(shí), 檢驗(yàn)勢(shì)接近于1, 平均運(yùn)行長(zhǎng)度為11.5, 誤報(bào)率為0; 但當(dāng)變點(diǎn)的位置離監(jiān)測(cè)起始時(shí)刻較遠(yuǎn)時(shí), 檢驗(yàn)勢(shì)迅速下降, 誤報(bào)率迅速增加, 此時(shí)監(jiān)測(cè)到的變點(diǎn)大都發(fā)生在m+k*之前, 所以其平均運(yùn)行長(zhǎng)度為負(fù)值.

總的來(lái)說(shuō), 統(tǒng)計(jì)量Stat在兩種邊界函數(shù)下的平均運(yùn)行長(zhǎng)度較長(zhǎng), 但誤報(bào)率很小, 大多數(shù)情況下都為0.CUSUM方法在變點(diǎn)發(fā)生在監(jiān)測(cè)起始時(shí)刻時(shí), ARL較短, 誤報(bào)率較低.但當(dāng)變點(diǎn)發(fā)生在離起始時(shí)刻較遠(yuǎn)時(shí), 誤報(bào)率較高.

5.實(shí)例分析

圖1是一組新生兒睡眠狀態(tài)的觀測(cè)序列, 每隔30秒觀測(cè)一次, 共1000個(gè)觀測(cè)值.每個(gè)時(shí)刻的睡眠狀態(tài)一般可分為四類: 安靜睡眠(Quiet sleep), 不確定睡眠(Indeterminate sleep), 主動(dòng)睡眠(Active sleep), 醒(Awake), 分別記為1, 2, 3, 4.根據(jù)新生兒的睡眠規(guī)律, 四個(gè)睡眠狀態(tài)有一定的次序, 即“4”＜“1”＜“2”＜“3”, 則觀測(cè)序列可看作有序分類時(shí)間序列.文[20]利用累積logistic回歸模型對(duì)其進(jìn)行建模

其中α1=-14.722,α2=-10.389,α3=-4.078,β1= 18.663,β2= 12.173,β3= 7.566.采用文[13]的回顧性變點(diǎn)檢驗(yàn)(Retrospective change-point detection)方法檢測(cè)上述模型是否存在結(jié)構(gòu)變點(diǎn), 結(jié)果發(fā)現(xiàn)α2在第596個(gè)觀測(cè)值處存在變化.在變點(diǎn)前后對(duì)模型參數(shù)重新進(jìn)行估計(jì), 得到調(diào)整模型.對(duì)比原模型和調(diào)整后模型的AIC值, 發(fā)現(xiàn)調(diào)整后的模型AIC值減小, 即調(diào)整后的模型更加合理.用本文所提出的變點(diǎn)序貫檢驗(yàn)監(jiān)測(cè)變點(diǎn), 以前200個(gè)數(shù)據(jù)為歷史數(shù)據(jù), 從第201個(gè)數(shù)據(jù)開始監(jiān)測(cè), 發(fā)現(xiàn)統(tǒng)計(jì)量在第671個(gè)數(shù)據(jù)處停止監(jiān)測(cè), 監(jiān)測(cè)延遲為175.雖然有一定的監(jiān)測(cè)延遲, 但也能有效監(jiān)測(cè)到變點(diǎn).

圖1 新生兒睡眠狀態(tài)的觀測(cè)序列

6.結(jié)論

累積Logistic回歸模型是分類時(shí)間序列建模的一類重要模型.針對(duì)累積Logistic回歸模型結(jié)構(gòu)變點(diǎn)的序貫檢驗(yàn)問(wèn)題, 本文基于偏似然得分過(guò)程構(gòu)造統(tǒng)計(jì)量, 在線監(jiān)測(cè)模型的結(jié)構(gòu)是否存在變點(diǎn).在原假設(shè)下推導(dǎo)出統(tǒng)計(jì)量的漸近分布, 備擇假設(shè)下證明其一致性.模擬試驗(yàn)采用了兩種形式的邊界函數(shù),結(jié)果表明原假設(shè)下統(tǒng)計(jì)量的檢驗(yàn)水平大部分都接近于顯著性水平0.05, 少部分存在一些扭曲.在備擇假設(shè)下雖然監(jiān)測(cè)到變點(diǎn)的平均運(yùn)行長(zhǎng)度較長(zhǎng), 但監(jiān)測(cè)到變點(diǎn)的誤報(bào)率較低, 且檢驗(yàn)勢(shì)也較高.CUSUM變點(diǎn)檢測(cè)方法在變點(diǎn)發(fā)生在監(jiān)測(cè)起始時(shí)刻時(shí), ARL較短, 誤報(bào)率較低.但當(dāng)變點(diǎn)發(fā)生在離起始時(shí)刻較遠(yuǎn)時(shí), 誤報(bào)率較高.最后通過(guò)一組新生兒睡眠狀態(tài)的觀測(cè)數(shù)據(jù)說(shuō)明本文方法的有效性.