一種基于Newton迭代法的累積Logistic回歸模型參數(shù)估計(jì)

印明昂， 王鈺爍， 孫志禮， 郭 兵

(東北大學(xué) 機(jī)械工程與自動(dòng)化學(xué)院， 遼寧 沈陽(yáng) 110819)

累積Logistic回歸又稱(chēng)為有序多分類(lèi)Logistic回歸、次序Logistic回歸[1]，其具有模型變量少、分類(lèi)效率高的特點(diǎn).近年來(lái)，該方法在醫(yī)藥、金融及制造等領(lǐng)域得到了廣泛的應(yīng)用.文獻(xiàn)[2]選取2015年3月至2017年4月接受培非格司亭(Pegfilgrastim)聯(lián)合化療的166例癌癥案例，采用多變量有序Logistic回歸分析，確定符合Pegfilgrastim預(yù)防以維持相對(duì)劑量強(qiáng)度的預(yù)測(cè)因素.文獻(xiàn)[3]基于發(fā)動(dòng)機(jī)健康變化的Logistic回歸模型，實(shí)現(xiàn)了貝葉斯?fàn)顟B(tài)估計(jì)的健康狀態(tài)序列更新，進(jìn)而預(yù)測(cè)引擎未來(lái)的健康變化.文獻(xiàn)[4]通過(guò)審查2003至2011年外來(lái)資本流入的綜合影響，檢驗(yàn)了次序Logistic回歸模型應(yīng)用于基于柯布-道格拉斯生產(chǎn)函數(shù)的實(shí)證經(jīng)濟(jì)增長(zhǎng)研究的可行性.

在應(yīng)用領(lǐng)域日益擴(kuò)展，數(shù)據(jù)數(shù)量和維度不斷增大的現(xiàn)狀下，目前，累積Logistic回歸模型的參數(shù)估計(jì)通常采用統(tǒng)計(jì)學(xué)軟件SPSS或SAS的自帶功能實(shí)現(xiàn).本文基于Newton迭代法提出一種累積Logistic回歸模型的參數(shù)估計(jì)方法，對(duì)其中比較敏感的迭代初值選取及迭代過(guò)程Hessian矩陣奇異問(wèn)題進(jìn)行了自適應(yīng)控制，最后利用美國(guó)凱斯西儲(chǔ)大學(xué)軸承數(shù)據(jù)庫(kù)(CWRU)[5]數(shù)據(jù)對(duì)該算法模型與SPSS模型進(jìn)行了對(duì)比驗(yàn)證.

1 參數(shù)估計(jì)

1.1 迭代格式的構(gòu)建

設(shè)(x,y)為有序J分類(lèi)樣本，其中輸入觀測(cè)值x為nM矩陣，n代表容量，M代表維度；有序輸出觀測(cè)值y為n1列向量，其元素yi,i∈{1,2,,n},yi可取值為1,2,,J.則第i個(gè)輸入xi=[xi1,xi2,,xiM]的前j(j∈{1,2,,J})類(lèi)累積Logistic表達(dá)式為[6]

(1)

其中：p(yi≤j|xi)為累積概率(簡(jiǎn)記為pij)，且有p(yi≤J|xi)=1；αj為截距，對(duì)于累積Logistic模型共有J-1個(gè)；βm為回歸系數(shù)，所有分類(lèi)的系數(shù)相同.由此可得

(2)

則該輸入屬于第j類(lèi)的概率為

(3)

最后由概率最大者確定其所屬類(lèi)別.不妨設(shè)Y為n×J矩陣，其元素yij為

(4)

則它的對(duì)數(shù)似然函數(shù)可表示為

(5)

由Newton迭代公式得到迭代格式：

bk+1=bk-H(bk)-1f(bk) .

(6)

(7)

(8)

H(bk)元素如下

pij-1(1-pij-1)]；

(9)

(10)

(11)

其中，r,j=1,2,,J-1;s,m=1,2,,M.

1.2 迭代初值的選取

Newton迭代法雖然具有平方收斂速度，但它對(duì)迭代初值非常敏感，如果選取不當(dāng)往往會(huì)導(dǎo)致迭代發(fā)散.下面通過(guò)數(shù)學(xué)推導(dǎo)，給出一種初值的選取方法.

首先將數(shù)據(jù)歸一化，使所有輸入數(shù)據(jù)統(tǒng)一映射到[0,1]區(qū)間上[7],即

(12)

其中：xim表示輸入i第m列的數(shù)據(jù)；xm為所有第m列數(shù)據(jù)組成的向量，i=1,2,,n，m=1,2,,M；min(·)，max(·)分別表示向量中的最小、最大值.

(13)

(14)

-fd1≤β0i≤fd1,i=1,2,,M.

(15)

由式(14)可推出：

(16)

其中，c=ln(1/p1-1) -ln(1/p2-1).進(jìn)而有

α10>>α(J-2)0>α(J-1)0.

(17)

式(14)～式(17)給出了b0的大致范圍.為進(jìn)一步壓縮區(qū)間，提高選取成功率，將b0分為三部分：b0=[α0T,β01T,β02T]T，其中α0=[α10,,α(J-1)0]T，β01=[β10,,βk0]T，β02=[β(k+1)0,,βM0]T，k[1,J].將式(14)的第一和最后一式進(jìn)行縮放：

(18)

(19)

其中fd2=M·fd1.設(shè)α0，β01已知，代入式(14)整理為關(guān)于β02的不等式組：

(21)

算法描述見(jiàn)算法1.

算法1 迭代初值的選取

輸入x，fd1，k.

2) 區(qū)間[-fd1,fd1]內(nèi)隨機(jī)生成β01;

3) 以式(19)計(jì)算區(qū)間，在該區(qū)間內(nèi)隨機(jī)生成α0;

4) 將α0降序處理，判斷式(16)是否成立.若是，轉(zhuǎn)步驟5)；若否，轉(zhuǎn)步驟2)；

5) 以式(15)，式(20)為約束，判斷式(21)是否存在最優(yōu)解.若是，得到β02，結(jié)束；若否，轉(zhuǎn)步驟2).

輸出b0=[α0T,β01T,β02T]T.

通過(guò)此算法得到的b0即可作為Newton法的迭代初值.

1.3 迭代過(guò)程的控制

雖然得到了較好的迭代初值，但是在迭代過(guò)程中可能會(huì)出現(xiàn)如下情況：

1)αk+1元素之間大小關(guān)系錯(cuò)亂；

2)βk+1與βk之間差距過(guò)大.

這些問(wèn)題會(huì)導(dǎo)致Hessian矩陣奇異，迭代失敗.通過(guò)分析，這些情況通常是由步長(zhǎng)Δk=H-1(bk)·f(bk) =[Δαk; Δβk]過(guò)大造成.下面給出一種自適應(yīng)控制方法.

當(dāng)出現(xiàn)情況1)時(shí)，利用下式減小步長(zhǎng)：

bk+1=bk-w·Δk.

(22)

其中，w為調(diào)整系數(shù)，初值為1，一旦常系數(shù)在k+1次迭代不滿(mǎn)足次序條件時(shí)便將w賦值為w/2，重新計(jì)算bk+1.

為防止情況2)的出現(xiàn)，需重點(diǎn)控制Δβk+1元素之間的差距.設(shè)置閾值td，令u=max(|Δβk+1|).若u大于閾值td，則需要對(duì)其進(jìn)行約束：

(23)

其中:放大系數(shù)fd3推薦為1.2～1.5；閾值td為1～2.

算法描述見(jiàn)算法2.

算法2 迭代過(guò)程的控制

輸入b0，fd3，td，收斂閾值ε.

1) 初始化參數(shù)k=0，w=1;

2) 利用迭代格式(6)計(jì)算bk+1;

3) 判斷αk+1元素間大小關(guān)系是否正確.若是，轉(zhuǎn)步驟4)；若否，w=w/2，通過(guò)式(22)調(diào)整，轉(zhuǎn)步驟3);

4) 判斷u=max(|Δβk+1|)是否在閾值td之內(nèi).若是，轉(zhuǎn)步驟5)；若否，通過(guò)式(23)調(diào)整，轉(zhuǎn)步驟3);

輸出bk+1.

算法2中的“=”為賦值符號(hào).由于算法1已保證了α0的順序，因此只要w調(diào)整適當(dāng)，總會(huì)保持αk的大小關(guān)系不變.

2 滾動(dòng)軸承健康狀態(tài)評(píng)估實(shí)例

文獻(xiàn)[9]通過(guò)非線性時(shí)間序列Logistic模型，驗(yàn)證了檢測(cè)滾動(dòng)軸承非線性突變信號(hào)的敏感性和有效性，因此本文也以軸承健康狀態(tài)評(píng)估為例.選取美國(guó)凱斯西儲(chǔ)大學(xué)軸承數(shù)據(jù)庫(kù)[5]中采樣頻率為12 kHz，負(fù)載為0，正常狀態(tài)和滾動(dòng)體故障(通過(guò)電火花加工設(shè)置故障，火花點(diǎn)直徑為0.177 8，0.355 6，0.533 4 mm，分別模擬輕微磨損、一般磨損及嚴(yán)重磨損)的振動(dòng)信號(hào)數(shù)據(jù)為原始數(shù)據(jù)，其中響應(yīng)類(lèi)別1，2，3，4分別表示正常狀態(tài)、輕微磨損、一般磨損和嚴(yán)重磨損，模型建立過(guò)程如下.

首先，根據(jù)文獻(xiàn)[10]，提取振動(dòng)信號(hào)中的30個(gè)特征，得到相應(yīng)特征值，再利用逐步模型選擇法，以Wald統(tǒng)計(jì)量為指標(biāo)從中逐步篩選出8個(gè)主要特征：FC，F(xiàn)2，F(xiàn)5，F(xiàn)6，F(xiàn)7，F(xiàn)8，F(xiàn)9，F(xiàn)12；然后，利用文獻(xiàn)[11]所述剔除離群點(diǎn)方法，從所有606組數(shù)據(jù)中篩選出有效數(shù)據(jù)597組；最后，通過(guò)本文方法和SPSS軟件分別建立累積Logistic模型，所得回歸系數(shù)如表1所示.

表1 不同模型的回歸系數(shù)

根據(jù)回歸結(jié)果，兩種方法所得各系數(shù)在組內(nèi)所占權(quán)重大致相同，且各自變量系數(shù)的正負(fù)符號(hào)相同，說(shuō)明兩種方法對(duì)于當(dāng)前數(shù)據(jù)具有相同的解釋性.模型評(píng)價(jià)指標(biāo)數(shù)值結(jié)果見(jiàn)表2.

表2 不同模型的評(píng)價(jià)指標(biāo)

表2中除赤池信息量AIC外，其余指標(biāo)的數(shù)值越大越能證明模型的優(yōu)勢(shì).可以看出，本文方法在各項(xiàng)指標(biāo)上都更加突出，證明了該方法的有效性.為進(jìn)一步驗(yàn)證穩(wěn)健性，分別利用兩種方法進(jìn)行200次Booststrap隨機(jī)試驗(yàn).每次試驗(yàn)在正常、輕微磨損、一般磨損及嚴(yán)重磨損這4種類(lèi)別中，分別隨機(jī)選取25組共100組數(shù)據(jù)作為驗(yàn)證樣本，其余作為訓(xùn)練樣本.將試驗(yàn)結(jié)果制成柱狀分布圖.兩方法的訓(xùn)練準(zhǔn)確率及驗(yàn)證準(zhǔn)確率如圖1～圖4所示.

由圖1～圖4可知，兩種方法訓(xùn)練準(zhǔn)確率和驗(yàn)證準(zhǔn)確率的分布都大致符合正態(tài)分布；通過(guò)對(duì)比可以直觀看出，本文方法無(wú)論訓(xùn)練準(zhǔn)確率還是驗(yàn)證準(zhǔn)確率，都處于更高水平，表明該方法具有較強(qiáng)的穩(wěn)健性.

3 結(jié) 論

1) 通過(guò)自適應(yīng)地選取迭代初值和控制迭代過(guò)程，避免了Hessian矩陣奇異的情況，使Newton迭代法能夠有效地進(jìn)行累積Logistic模型的參數(shù)估計(jì).

2) 與SPSS累積Logistic回歸模型數(shù)值結(jié)果對(duì)比，本文方法在各項(xiàng)模型評(píng)價(jià)指標(biāo)上具有更高的水平，驗(yàn)證了該方法的有效性.

3) 基于200次的Booststrap隨機(jī)試驗(yàn)對(duì)比，表明本文方法對(duì)于數(shù)據(jù)的依賴(lài)性較小，具有較強(qiáng)的穩(wěn)健性.