鄧文麗,朱瑩瑩
(江西師范大學(xué) 數(shù)學(xué)與信息科學(xué)學(xué)院,南昌 330027)
約束條件下的統(tǒng)計(jì)推斷已經(jīng)成為統(tǒng)計(jì)分析的一個(gè)重要領(lǐng)域,保序回歸是約束條件下的統(tǒng)計(jì)推斷的一類典型問題。保序回歸的研究從20世紀(jì)中葉就已經(jīng)開始,到20世紀(jì)60年代,保序回歸得到了廣泛的重視。Bartholoremew,Barlow等對(duì)此進(jìn)行了研究,并且于1972年聯(lián)合出版了保序回歸的第一本專著 《The theory and application of Isotonic Regression》,使得保序回歸問題得到了進(jìn)一步的重視和發(fā)展,成為熱門話題。
在回歸分析問題中,假定自變量為 z=(z1,z2,…,zn),因變量為 y=(y1,y2,…,yn)。 在很多情況下,研究者并沒有關(guān)于回歸函數(shù)的很多信息,而只是一些關(guān)于變量的約束條件,比較典型的一類就是,假設(shè)z的元素是經(jīng)過排序的,z1≤z2≤…≤zn,yi會(huì)隨著zi遞增,這樣一種情形被稱為保序回歸;yi會(huì)隨著zi遞減,被稱為逆序回歸。兩種情況合在一起,被稱為單調(diào)回歸。詳細(xì)介紹可參見文獻(xiàn) de Leeuw(2005)。
令 x=(x1,x2,…,xn)表示 y 所對(duì)應(yīng)的回歸擬合值,保序回歸可以表述為:
尋找 x*=(x1*,z2*,…,zn*)∈G,G={x∈Rn|x1≤x2≤…≤xn}使得
其中 ω=(ω1,ω2,…,ωn),ωi≥0 (i=1,2,…,k)是給定的權(quán)數(shù),ω1+ω2+…+ωn=1。這時(shí)x*被稱為y的保序回歸。
保序回歸問題的求解過程中,有一個(gè)最基本的定理是:如果yi>yi+1,那么它們對(duì)應(yīng)的擬合值滿足x^i=x^i+1。因此,保序回歸的求解歸結(jié)為找到指標(biāo)集{1,2,…,n}的一個(gè)合理劃分,B1,B2,…Bk,k≤n,使得 B1+B2+…+Bk={1,2,…,n},對(duì)劃分中的任意一塊Bj而言,任給i∈Bj,xi都是相等的,在滿足這樣條件的劃分中,尋找(1)的解。 Ayer et al.(1995)提出的 PAVA(Pool Adjecent Violators Algorithm)算法使問題(1)從理論上得到了很好的解決。PAVA算法求得的“最優(yōu)”劃分B1,B2,…Bk,k≤n,對(duì)劃分中的任意一塊 Bj而言,任給 i∈Jj,。
本文將經(jīng)典的保序回歸問題拓展到絕對(duì)損失下,去解決下面的最優(yōu)化問題:
尋找 x*∈G,G={x∈Rn|x1≤x2≤…≤xn}使得
其中 ω=(ω1,ω2,…,ωn),ωi≥0 (i=1,2,…,k)是給定的權(quán)數(shù),ω1+ω2+…+ωn=1。
本文將根據(jù)PAVA思想過程給出一種迭代算法,通過迭代過程分析得出最優(yōu)化問題(2)的解,并給予證明。
J=(B1,B2,…,Bk),B1,B2,…,Bk,k≤n 是指標(biāo)集{1,2,…,n}的一個(gè)劃分。 若 B={p,p+1,…,q},1≤p≤q≤n,J中包含 p-1的元素,用B-表示;包含q+1的元素用B+表示。若p=1時(shí),B-=φ,若 q=n 時(shí),則 B+=φ。 用[AM(B)L,AM(B)U]記指標(biāo)集 B所對(duì)應(yīng)y的元素的加權(quán)中位數(shù)所在集合,有
類似PAVA算法的思想,下面針對(duì)本文討論的最優(yōu)化問題(2)給出一種迭代算法。
G 表示可行域{x1,x2,…,xn|x1≤x2≤…≤xn},如果(y1,y2,…,yn)∈G,則=yi,i=1,2,…,n 否則,令 J={{1},{2},…,{n}},B={1},進(jìn)入迭代算法。
(1)若 B+=φ,則迭代結(jié)束,xi*∈[AM(B)L,AM(B)U],i∈B,B∈J;否則將 AM(B)L,AM(B)U和 AM(B+)L,AM(B+)U進(jìn)行比較。
若 AM(B)U≤AM(B+)L,B=B,回到初始;
若 AM(B)U>AM(B+)L,但 AM(B)L≤AM(B+)U,即
若 AM(B)L>AM(B+)U進(jìn)入第(2)步。
(2)J=J{B,B+}∪(B∪B+),令 B=B∪B+,計(jì)算AM(B)L=AM(B)U。
(3)若 B-=φ,進(jìn)入第(1)步;否則進(jìn)行下面的判斷:
若 AM(B-)U≤AM(B)L,回到第(1)步;
若 AM(B-)U>AM(B)L,但 AM(B-)L≤AM(B)U,即
進(jìn)入第(1)步。
AM(B-)L>AM(B)U,進(jìn)入第(4)步;
(4)J=J{B,B-}∪(B∪B-),B=B∪B+,進(jìn)入第(3)步。
迭代的結(jié)果,指標(biāo)集{1,2,…,n}被劃分成了 B1,B2,…,Bk,k≤n 是指標(biāo)集的一個(gè)劃分。 B1+B2+…+Bk={1,2,…,n},而且最優(yōu)解屬于集合
下面證明(5)所表示的就是最優(yōu)化問題(2)的解集。
定理 1 假定指標(biāo)集 B=(p,p+1,…,q),?i∈B,yi∈R,yi所對(duì)應(yīng)的權(quán)數(shù)記為 wi,yp≥yp+1≥…≥yq。 數(shù)列{yi,i∈B}的加權(quán)中位數(shù)所在的集合表示為[AM(B)L,AM(B)U],[AM(B)L,AM(B)U]如(3)、(4)所示,那么
①使得
定理 2 若 B=(p,p+1, …,q),1≤p≤q≤n 任給 p≤i≤+1≥yi+2≥…≥yq,記
①若 T≠φ,選取
可使得
證明:①式顯然成立,下面證②。
根據(jù)定理1中的結(jié)論①,又因?yàn)閥p≥yp+1≥…≥yi,yi+1≥yi+2≥…≥yq, 那么使得成立的一定滿足使得成立的一定滿足…=。 所以使得(7)式成立的滿足
因?yàn)?T=φ,即 AM(Ui(B))U<AM(Li(B))L,顯然使得(7)式成立的{,i∈B}還應(yīng)滿足
定理 3 若 B={p,p+1,…q},1≤p≤q≤n,任給 p≤i≤q,Li(B){p,p+1,…,i}。 記記
①若T≠φ,選取
可使得
例:已知 y1=4,y2=5,y3=1,y4=6,y5=8,y6=7 并且 wi=1/6,i=1,…,6。 求 x* 使在滿足 x1≤x2≤…≤x6的限制下達(dá)到最小。
解:y1=4<y2=5,y2=5>y3=1,所以 B1={1},B2={2,3}
取 AM(B1)L=AM(B1)U=4,AM(B2)L=AM(B2)U=5,由此重新記 AM(B1)L=AM(B1)U=4,AM(B2)L=4,AM(B2)U=5;AM(B2)U<y4<y5;而 y5>y6,所以 B3={4},AM(B3)L=AM(B4)U=6;B4={5,6},AM(B4)L=AM(B4)U=8,顯然 AM(B4)L>AM(B3)U=6。 所以選取,J={{1},{2,3},{4},{5,6}}[4,5],∈[7,8],
用上述迭代計(jì)算該模型的解的次數(shù)為O(n),當(dāng)n不是很大時(shí),其運(yùn)算并不復(fù)雜,當(dāng)n很大時(shí),利用計(jì)算機(jī)輔助仍然可以比較快得出結(jié)果。因此,加權(quán)最小絕對(duì)偏差在保序限制下迭代計(jì)算仍是比較實(shí)用的,有奇異值出現(xiàn)時(shí),其所得解也比加權(quán)最小平方損失下的解更穩(wěn)健。
[1]Michael J BEST,Nilotpal Chakravarti.Active Set Algorithms for Isotonic Regression;A Unifying Framework[J].Mathematical Programming,1990,(47).
[2]De Leeuw,Jan,Hornik,Kurt,Mair,Patrick.Isotone Optimization in R:Pool-Adjacent-Violators Algorithm (PAVA) and Active Set Methods[M].UC Los Angeles:Department of Statistics,2009.
[3]De Leeuw.Monotonic Regression,Encyclopedia of Statistic in Behavioral Science[M].New York:Wiely,2005.
[4]史寧中.保序回歸與最大似然估計(jì)[J].應(yīng)用概率統(tǒng)計(jì),1993,(19).
[5]刑務(wù)強(qiáng).保序回歸的研究及應(yīng)用[D].西北工業(yè)大學(xué),2002.
[6]趙選民,刑務(wù)強(qiáng).凸規(guī)劃下的保序回歸[J].數(shù)理統(tǒng)計(jì)與管理,2003,(22).