楊豐凱，袁海靜

0 引言

變點(diǎn)問(wèn)題是指在一隨機(jī)序列中存在某一時(shí)刻，使得該時(shí)刻兩側(cè)的序列服從不同的分布。

從上世紀(jì)50年代開(kāi)始，估計(jì)一隨機(jī)序列中變點(diǎn)的位置成為統(tǒng)計(jì)學(xué)中的研究熱點(diǎn)之一，Chen等[1]詳細(xì)介紹了各種變點(diǎn)模型及其在遺傳學(xué)、醫(yī)藥以及金融領(lǐng)域的應(yīng)用。其中一類(lèi)問(wèn)題是研究如何有效地估計(jì)線(xiàn)性回歸模型中回歸系數(shù)的變點(diǎn)位置，該類(lèi)問(wèn)題可描述為：對(duì)于序列yi，i=1，…n， 存在位置 r，p≤r≤n-p，使得：

且 εi，i=1，…，n 相 互 獨(dú) 立 。 其 中 xi=(1，xi1，…xi，p-1)T,β1，β2為不同的 p 維回歸系數(shù)。本文的目標(biāo)是估計(jì)變點(diǎn)位置r。對(duì)于該類(lèi)問(wèn)題，文獻(xiàn)中的研究方法多是基于似然的方法和貝葉斯方法。其中，Quandt[2,3]發(fā)展了基于最大似然估計(jì)和似然比檢驗(yàn)的回歸系數(shù)變點(diǎn)檢驗(yàn)和估計(jì)方法，Kim等[4,5]研究了回歸系數(shù)變點(diǎn)似然比檢驗(yàn)的漸進(jìn)性質(zhì)。 Ferreira[6]，Chin Choy 和 Broemeling[7]，Holbert[8]則詳細(xì)討論了回歸系數(shù)變點(diǎn)估計(jì)的貝葉斯方法。Chen等[1,9]則從信息論的角度，提出了基于Schwarz信息準(zhǔn)則的變點(diǎn)估計(jì)方法?；隈R氏鏈蒙特卡洛（MCMC）的Gibbs抽樣，由于其靈活性和易實(shí)施性，是一種有效的貝葉斯變點(diǎn)估計(jì)方法，但Gibbs抽樣是一種迭代抽樣算法，所抽取的馬氏鏈?zhǔn)欠袷諗康胶篁?yàn)分布很難判斷，并且所抽取的樣本也很難保證是獨(dú)立的。Tian等[10]發(fā)展了一種基于逆貝葉斯公式的非迭代抽樣算法，稱(chēng)為IBF抽樣，該算法能夠直接從離散的后驗(yàn)分布中抽取獨(dú)立同分布的樣本，然后依據(jù)該樣本對(duì)相關(guān)參數(shù)做統(tǒng)計(jì)推斷，從而巧妙地避開(kāi)了Gibbs抽樣的不足之處。Tian等[11]將IBF抽樣算法應(yīng)用到變點(diǎn)問(wèn)題的研究中，并討論了泊松變點(diǎn)在醫(yī)學(xué)數(shù)據(jù)中的應(yīng)用。Yang等[12]則研究了正態(tài)分布均值變點(diǎn)的IBF抽樣算法。本文基于Yang等[12]的研究，將IBF算法應(yīng)用到模型（1）中，討論估計(jì)回歸系數(shù)變點(diǎn)位置的非迭代抽樣算法。分別在弱先驗(yàn)信息和共軛先驗(yàn)信息下給出相應(yīng)變點(diǎn)位置的精確后驗(yàn)分布，用IBF獲得獨(dú)立同分布的樣本，并作相應(yīng)的統(tǒng)計(jì)推斷，大量模擬結(jié)果顯示，該算法能夠有效地估計(jì)變點(diǎn)位置，并且算法的運(yùn)行速度比迭代的Gibbs抽樣要快很多。

1 無(wú)信息先驗(yàn)下正態(tài)回歸均值變點(diǎn)估計(jì)的IBF算法

本文在無(wú)信息先驗(yàn)分布下討論回歸系數(shù)變點(diǎn)估計(jì)的非迭代抽樣算法。記 y=(y1，…，yn)T，模型（1）的似然函數(shù)為：

且 β1，β2，σ2，r 相互獨(dú)立。 其中 U(p，…，n-p)表示 p～(n-p)上的離散均勻分布，“∝”為正比符號(hào)。

1.1 無(wú)信息先驗(yàn)下的條件分布

參數(shù) (β1，β2，σ2，r)的聯(lián)合后驗(yàn)分布：

取如下無(wú)信息先驗(yàn)分布：

從而變點(diǎn)位置r的條件后驗(yàn)分布為：

而：

上式右端為逆伽馬分布的密度的核，故有：

其中 IG(a，b)表示參數(shù)為a，b的逆伽馬分布。

1.2 IBF算法

根據(jù)Tian等[10,11]所提出的離散缺失數(shù)據(jù)分析的IBF算法的思想，有：

IBF抽樣算法分以下幾步：

第一步：計(jì)算

則變點(diǎn)位置r的精確后驗(yàn)分布為：

P(r=i|y)=λi，i=p，…，n-p

第二步：以概率 λ=(λp，…，λn-p)從集合 S={p，…，布π(r|y)的獨(dú)立同分布的樣本，可以根據(jù)該樣本對(duì)r作統(tǒng)計(jì)推斷。

第三步：對(duì) l=1，…，L,根據(jù)式（11）、式（6）和式（9），分別產(chǎn)生：

2 共軛先驗(yàn)下正態(tài)回歸均值變點(diǎn)估計(jì)的IBF算法

在共軛先驗(yàn)下研究正態(tài)回歸均值變點(diǎn)估計(jì)的IBF算法，取如下先驗(yàn)分布：

其中μ1、μ2均為 p維常向量，D1、D2均為 p階正定陣，IG(a2，b2)表示參數(shù)為a2和b2的逆伽瑪分布。在此假設(shè) β1，β2關(guān)于 σ2條件獨(dú)立,且 (β1，β2，σ2)與 r 相互獨(dú)立。

2.1 共軛先驗(yàn)下的條件分布

參數(shù) (β1，β2，σ2，r)的聯(lián)合后驗(yàn)分布：

從而變點(diǎn)位置r的條件后驗(yàn)分布為：

該分布與式（2）相同，與先驗(yàn)分布的選取無(wú)關(guān)。而：

分別計(jì)算上式等號(hào)右端的三個(gè)條件分布。易見(jiàn)：

2.2 IBF算法

類(lèi)似無(wú)信息先驗(yàn)下的IBF算法，有：

步驟3：對(duì)l=1，…，L,根據(jù)式（19）、式（15）和式（17），分別產(chǎn)生：

其中：

3 模擬

通過(guò)模擬，研究用IBF算法估計(jì)回歸系數(shù)變點(diǎn)位置的效果，并與迭代的Gibbs抽樣算法作比較。考慮如下單變點(diǎn)一元線(xiàn)性回歸模型：

假設(shè)εi，i=1，…，n相互獨(dú)立。在此取n=200，β10=選取變點(diǎn)位置r=40，60，80，100，120 ， 140，160。重復(fù)試驗(yàn)200次，在每次實(shí)驗(yàn)中分別實(shí)施IBF算法和Gibbs抽樣算法對(duì)變點(diǎn)位置進(jìn)行估計(jì)。在第i次試驗(yàn)中，記：

分別取前一半數(shù)據(jù)和后一半數(shù)據(jù)所得最小二乘估計(jì)作為初值，即有：

取σ2(0)=1。在Gibbs抽樣中，取r(0)=100。采取基于樣本的共軛先驗(yàn)分布，其中：

分別模擬產(chǎn)生變點(diǎn)位置r的L=6000個(gè)IBF樣本和L=6000個(gè)Gibbs樣本，并舍去前3000個(gè)作為burn-in樣本，方誤差的平方根（RMS）來(lái)衡量估計(jì)的精度，即有：

模擬結(jié)果如表1所示。

表1 對(duì)不同的變點(diǎn)位置，兩種算法的結(jié)果比較

表1顯示不論真實(shí)變點(diǎn)在何位置，兩種算法所得估計(jì)r?都很接近真實(shí)變點(diǎn)位置，并且估計(jì)的精度RMS都很小，說(shuō)明兩種算法對(duì)真實(shí)的變點(diǎn)位置不敏感，都能夠有效地估計(jì)變點(diǎn)位置，并且兩種算法所得估計(jì)r?以及估計(jì)的精度RMS相差微乎其微。

比較對(duì)于產(chǎn)生相同的樣本量的兩種算法的計(jì)算時(shí)間。以r=80為例，本文取J分別為200,400,800,1000,3000,6000,10000,20000。共做200次重復(fù)試驗(yàn)，這200次試驗(yàn)的平均運(yùn)行時(shí)間time(s)及r的估計(jì)r?和均方誤差的平方根（RMS）見(jiàn)表2。

表2 對(duì)不同的樣本量?jī)煞N算法的結(jié)果比較

從表2可以看出，對(duì)不同的樣本量，兩種算法所得變點(diǎn)位置的估計(jì)值都很接近真值，所得估計(jì)的RMS都很小，說(shuō)明這兩種算法都很有效。但兩種算法所產(chǎn)生相同樣本量的運(yùn)行時(shí)間卻相差很大。當(dāng)J從200增加到20000時(shí)，200次試驗(yàn)中IBF算法的平均運(yùn)行時(shí)間在0.0316秒到0.0325秒之間，相差只有0.0009秒，變化微乎其微；與之對(duì)比，Gibbs抽樣的平均運(yùn)行時(shí)間從0.965秒增大到99.399秒，后者是前者的大約100倍?？v向來(lái)比，產(chǎn)生J個(gè)樣本的時(shí)間比Gibbs/IBF，當(dāng)J=200時(shí)，該比值為30.54；當(dāng)J=20000時(shí)，該比值增大到3058.43。說(shuō)明要產(chǎn)生相同的樣本量J=20000，Gibbs抽樣的運(yùn)行時(shí)間是IBF算法的運(yùn)行時(shí)間的3058.43倍?？梢?jiàn)IBF抽樣比Gibbs抽樣快很多，特別是產(chǎn)生大樣本時(shí)。

4 結(jié)論

本文提出了估計(jì)線(xiàn)性回歸模型中回歸系數(shù)變點(diǎn)位置的非迭代抽樣算法（IBF），該算法能夠獲得變點(diǎn)位置的精確后驗(yàn)分布，進(jìn)而得到該后驗(yàn)分布的的獨(dú)立同分布的樣本,然后依據(jù)該樣本對(duì)變點(diǎn)位置做統(tǒng)計(jì)推斷。該算法巧妙地避開(kāi)了Gibbs抽樣等MCMC方法的收斂性診斷問(wèn)題，所獲樣本為簡(jiǎn)單隨機(jī)樣本，可直接用來(lái)進(jìn)行統(tǒng)計(jì)推斷。模擬顯示該算法能夠有效地估計(jì)未知變點(diǎn)位置，并且與迭代的Gibbs抽樣相比，該非迭代抽樣算法的運(yùn)行時(shí)間大大縮短。

[1]Chen J,Gupta A K.Parametric Statistical Change Point Analysis:With Applications to Genetics,Medicine,and Finance(2nd edition)[M].Boston:Birkhauser,2012.

[2]Quandt R E.The Estimation of The Parameters of a Linear Regression System Obeys Two Separate Regimes[J].Journal of the American Statistical Association,1958,(53).

[3]Quandt R E.Tests of the Hypothesis That a Linear Regression System Obeys Two Separate Regimes[J].Journal of the American Statistical Association,1960,(55).

[4]Kim H J,Siegmund D.The Likelihood Ratio Test for a Change-Point in Simple Linear Regression[J].Biometrika,1989,76(3).

[5]Kim H J.Tests for a Change-point in Linear Regression[J].IMS Lecture Notes-Monograph Series,1994,(23).

[6]Ferreira P E.A Bayesian Analysis of a Switching Regression Model:Known Number of Regimes[J].Journal of the American Statistical Association,1975,(70).

[7]Chin Choy J H,Broemeling L D.Some Bayesian Inferences for a Changing Linear Model[J].Technometrics,1980,(22).

[8]Holbert D.A Bayesian Analysis of a Switching Linear Model[J].Journal of Econometrics,1982,(19).

[9]Chen J.Testing for a Change Point in Linear Regression Models[J].Communications in Statistics-Theory and Methods,1998,(27).

[10]Tian G L,Tan M,Ng K W.An Exact Non-iterative Sampling Procedure for Discrete Missing Data Problems[J].Statistica Neerlandica,2007,61(2).

[11]Tian G L,Ng K W,Li K C et al.Non-iterative Sampling-based Bayesian Methods for Identifying Change Points in the Sequence of Cases of Hemolytic Uremic Syndrome[J].Computational Statistics and Data Analysis,2009,53(9).

[12]楊豐凱，袁海靜.正態(tài)均值變點(diǎn)識(shí)別的非迭代抽樣算法[J].統(tǒng)計(jì)與決策,2016，（8）.