一種基于廣義極值分布的非平衡數據分類算法

付俊杰 劉功申

(上海交通大學電子信息與電氣工程學院 上海 200240) (tianzhiyinyi@sjtu.edu.cn)

二分類是在不同領域的諸多應用中的一項非常常見的任務，如欺詐檢測、疾病預測、風險管理、文本分類等.在這些應用中，其中一個類別的數據數量遠超出另一個類別，使得傳統的機器學習方法難以得到理想的結果[1].除了直接預測數據對應的分類標簽，許多應用還可能關心這個預測的準確性有多少，也就是想知道數據屬于某個分類的概率[2].相對直接得到分類標簽，分類概率的預測往往能為人們的決策判斷提供額外的幫助.比如在醫療領域，醫生經常會通過聲學成像得到的圖像數據來判斷病人是否患有某種疾病，比起直接診斷病人是否患病，能說出病人有百分之幾的可能性患病將更有說服力[3].

非平衡數據分類的常用方法包括過采樣/欠采樣方法[4]、代價敏感學習[5]和集成學習，比如boosting方法[6].其中，采樣方法比較簡單且適應性較強，該方法的出發點簡單清晰，既然2個類別的數據量差異懸殊，就想辦法將它們拉近，然后用經典的分類算法解決問題；代價敏感學習與直接考慮數據數量的采樣方法不同，它考慮的是數據分錯類別時的代價，使分類器能更看重少類數據；集成學習先將原數據集劃分為多個子集，然后分別對每個子集訓練得到多個子分類器，最后將這些子分類器集成得到一個分類器，提高少類的分類準確性.然而到現在為止，關注分類概率預測的研究還較少.為了解決這個問題，一個自然的選擇便是廣義線性模型[7].眾所周知，廣義線性模型具有模型簡單、魯棒性強、泛用性廣等優點，在不同的鏈接函數下，廣義線性模型展現出了其良好的靈活性，更重要的是，該模型能直接給出分類的概率預測.然而，廣義線性模型中通常使用的鏈接函數對非平衡數據的適應性不佳[8]，因為大多數鏈接函數的曲線形狀是對稱的，從而導致它們對不同的錯分給予了相同的代價.因此，有研究表明，用非對稱的鏈接函數將更有利于非平衡數據分類[9-10]，而其中一種非對稱鏈接函數便是廣義極值(generalized extreme value, GEV)分布的累積分布函數.由于形狀參數ξ的存在， GEV鏈接函數能以不同的曲線斜率提供豐富多樣的曲線類型，可以適應現實生活中的各種不同平衡度的數據.

GEV分布廣泛應用于金融、氣象、計算機視覺等領域.近年來，將GEV分布作為鏈接函數運用在線性模型中的做法受到了一定的關注：Wang等人[9]將GEV鏈接函數運用到B2B電子支付系統的二分類問題中，得到了很好的結果；Calabrese等人[10]則在信貸領域的稀有事件分類中運用了GEV分布，通過最大化似然概率得到了模型的解；Agarwal等人[11]更是根據GEV分布的形式求得了對應的標準優化函數，使整個問題化為凸優化問題，在精確的分類概率預測問題上給出了一種解決方案.在這些研究中所用到的數據均是非平衡數據，而且均采用了GEV鏈接函數，這表明GEV分布在解決非平衡數據分類問題中具有一定潛力.

1 GEV線性回歸

1.1 廣義線性模型與GEV分布

令y∈{0,1}表示數據的分類標簽，x=(1,α1,α2,…,αd)T∈d+1表示輸入的數據，給定x的條件下y的期望為

E[y|x]=P(y=1|x)=g(βTx)，

(1)

其中,β是權重參數，g的逆函數被稱為鏈接函數.鏈接函數可以是任意的嚴格遞增函數，其定義域與值域符合g-1:[0,1]→，常用的有Logit鏈接、Probit鏈接和Cloglog鏈接等，前兩者是對稱鏈接函數，而Cloglog鏈接則是非對稱的.

GEV鏈接函數也是非對稱函數，與Cloglog不同的是，GEV鏈接函數具有可以調節非對稱程度的形狀參數.當均值為0、方差為1時，GEV分布的累積分布函數為

(2)

其中(·)+表示max(0,·);ξ是形狀參數，控制著函數曲線形狀以及定義域.GEV鏈接函數和其他的鏈接函數有2個主要區別：1)GEV鏈接函數具有可調的形狀參數，使得GEV鏈接函數可以通過改變曲線形狀來適應不同的數據，尤其是具有不同非平衡程度的數據；2)GEV鏈接函數族中的大多數曲線都是非對稱的，圖1展示了在不同的形狀參數下GEV分布的累積分布函數(CDF)和概率密度函數曲線(PDF).曲線的非對稱性意味著GEV鏈接函數對不同的分類錯誤會給予不同的錯分代價，而這正是處理非平衡數據分類問題的一個有效策略.

Fig. 1 CDF and PDF of GEV distribution with different values of ξ圖1 不同ξ值的GEV分布的累積分布函數和概率密度函數

1.2 GEV線性回歸

在廣義線性模型框架下，為了求得權重參數β，通常會以最大化似然概率為目標.假設有n個獨立同分布的數據樣本，那么廣義線性模型的目標優化函數即是：

(3)

但是當形狀參數ξ處在特定的區間時，上面的目標優化函數并非凸函數，也即無法求得全局最優解.為了解決這個問題，Agarwal等人[11]利用組合損失函數的方法求得了GEV鏈接函數對應的標準損失函數，使其配對形成凸函數優化.該方法在小數據集上有較好的表現，但由于其每輪運算中都需要計算海森矩陣的逆矩陣，該方法在大數據集上有計算性能問題.

為了在保證凸函數優化的基礎上進一步提升運算性能，本文運用了校準損失函數作為目標優化函數，該目標優化函數能以更通用的形式支持許多鏈接函數，使其能組合成為凸函數[12].校準損失函數為

l(β)=G(βTx)-yβTx，

(4)

其中，G為逆鏈接函數g的原函數.在這里g必須是單調遞增的以保證原函數的存在以及該函數的凸性.而原函數G不必直接求出，因為在算法迭代中需要計算的是梯度和海森矩陣，它們分別是：

(5)

(6)

因為g是單調遞增的，也即g′(βTx)≥0，容易證明海森矩陣是半正定的，因此最小化上面的校準損失函數可以轉換成凸優化問題.

如果函數g滿足利普希茨連續性，那么迭代算法的效率就能通過近似計算海森矩陣而得到提升.利普希茨連續性意味著存在一個常數L≥0，使得對于所有u，v∈，以下關系成立：

|g(u)-g(v)|≤L|u-v|，

(7)

式(7)表明函數g的斜率有最大值L.如果這個斜率最大值L能計算出來，那么我們可以簡單地將其替換到g′(βTx)的位置，使得牛頓法的迭代過程簡化為

βnew=βold-(LxxT)-1·(g(βTx)-y)x|β=βold.

(8)

此處用斜率最大值L替換精確的梯度，相當于在梯度更新方向正確的前提下提供了最小步長，而各維的實際更新步長仍受原始數據xxT的影響，保證了一定的參數更新效率.在替換過后，每輪的迭代計算中逆矩陣的運算與β無關，逆矩陣的計算可以放在循環開始前，這樣就能大大提高計算的效率.這種利普希茨優化既保留了牛頓法中的原始數據二階信息結構，還能像梯度下降法那樣有效率地更新，繼承了兩者的優點.

對于GEV鏈接函數，其符合利普希茨連續性，斜率最大值L為

L=exp(-(1+ξ))(1+ξ)1+ξ，

(9)

其中，形狀參數ξ≥-1，否則g′(βTx)會趨向無窮，有意義的L值將不存在.雖然這限制了ξ的取值，但有研究表明，ξ的取值在[-0.5,0.5]之間就能提供足夠的曲線形狀來使GEV鏈接函數具有高度的靈活性[9].

完整的算法如算法1所示.為避免過擬合，算法1中添加了L-2正則項.形狀參數ξ是一個重要的參數，它直接影響模型對非平衡數據的擬合程度，所以關于它的選取方法將放到第2節單獨介紹.β的初始值可以隨機設為接近0的值，或通過第2節介紹的方法設置.

算法1. GEV線性回歸算法.

輸入：XT=(x0,x1,…,xn)，xi∈k，Y=(y1,y2,…,yn)T，yi∈{0,1}，使用概率模型選優得到GEV形狀參數ξ∈[-1,+∞)，正則參數λ>0；

輸出：權重參數β.

初始化：

通過式(9)計算L.

步驟1. 計算近似海森矩陣的逆矩陣Σ-1=(L·XTXn+λI)-1，設置迭代下標t=0；

步驟2. 重復步驟2.1～2.4直到β收斂：

步驟2.3.βt+1=βt-Σ-1·zt；

步驟2.4.t=t+1.

2 形狀參數求取

為了能更好地發揮GEV鏈接函數的性能，形狀參數的選取尤為重要.交叉驗證可以作為一種備選方法，然而當需要選優的參數較多時，交叉驗證將變得非常耗時，因此這里使用概率模型對形狀參數選優.

2.1 最大化后驗概率

模型主要有2個參數待求解：權重參數β和形狀參數ξ.觀察到實際中形狀參數合適值在0附近，所以假設其服從標準正態分布，即P(ξ)=N(ξ|0,1).而權重參數各維數值不應過大，假設其各維度也服從標準正態分布，即P(βi)=N(βi|0,1)，并假設權重參數與形狀參數互相獨立.因此整個模型的后驗概率為

P(β,ξ|X,Y)∝P(Y|X,β,ξ)·P(β)·P(ξ)，

(10)

其中P(Y|X,β,ξ)是模型的似然概率，其單個yi在廣義線性模型框架下服從伯努利分布.對式(10)取對數后有：

(11)

Wang等人[9]在他們的研究中也用了類似的模型并通過最大化后驗概率(MAP)求得解，所不同的是他們還研究了其他的先驗概率分布，并且求解算法用的是馬爾可夫鏈蒙特卡羅方法.本文運用MAP方法是為了根據輸入數據自動選擇較優的形狀參數ξ.因此，模型中對權重參數β的求解不是必須的，權重參數是利用本文提出的GEV線性回歸算法來求解的.考慮到這一步驟的主要目標，本文對MAP的求解將采用迭代類算法，如牛頓法或梯度下降法，在每輪迭代中輪流更新權重參數和形狀參數.MAP目標優化函數分別對β和ξ求一階偏導數有：

(12)

(13)

其中,ai=1+ξηi.有了一階偏導數后就可以采用梯度下降法求取MAP目標優化函數的最大值，或繼續計算各自的二階偏導數利用牛頓法更新參數.得益于ξ的低維度，只需要很少的迭代次數就可以使形狀參數收斂到一個可接受的范圍.算法運行完后的副產物β參數可以看作是尋優過程的“中途解”，將之用作GEV線性回歸算法的初始值，可以進一步減少讓β參數收斂所需的迭代次數，改進回歸算法的效率.

2.2 貝葉斯分析

在MAP概率模型中，目標優化函數是后驗概率，我們假設GEV線性回歸算法中的形狀參數ξ和權重參數β都服從于一個確定的標準正態分布，這有時會變得不合適.貝葉斯分析中，我們依舊假設形狀參數ξ和權重參數β都服從正態分布，但是正態分布的2個關鍵參數——均值和方差，可以是未知的.為了方便計算，假設P(ξ|μξ)=N(ξ|μξ,1),P(β|μβ)=N(β|μβ,I),其中μξ和μβ分別是ξ的均值和β的均值，I是單位矩陣.通過證據近似法最大化邊緣似然概率P(Y|μξ,μβ)可以得到ξ和β.在計算中，ξ和β可以看作是隱變量，通過EM算法迭代求解.M步驟中需要最大化的目標函數Q為

Q(μξ,μβ|μξ,μβ)=

(14)

然而式(14)涉及對GEV分布的積分，直接計算很困難，所以這里通過Metropolis-Hastings采樣方法對此積分進行近似計算:

(15)

(16)

得到求和形式的函數Q后即可通過求導計算更新后的μξ和μβ,最優形狀參數ξ直接取μξ的均值即可.像MAP方法一樣，副產物μβ可以用于GEV線性回歸算法的初始化中.

3 實 驗

本節實驗主要分為3個部分：1)通過人工數據集驗證形狀參數求取方法的準確性；2)在真實數據集上橫向對比本文算法和其他線性回歸算法，在分類準確性和概率預測準確性上分別做出評價；3)對本文算法的運行效率進行實驗.

3.1 數據集與評價標準

人工數據集方面，數據維度為2維，采用互相獨立的標準正態分布隨機生成數據點xi，并固定權重參數β=(1,1)T，通過GEV逆鏈接函數生成分類概率值yi.本次實驗用的生成數據集使用了3個形狀參數值ξ，分別是-0.1，0.2，0.5；生成數據集的大小也有3種：100個、500個和1 000個，共9個生成數據集.真實數據集來源于Keel數據倉庫中的非平衡數據集[13].本次實驗根據數據所屬領域以及非平衡程度，從中挑出了17個真實數據集，它們均來自不同領域并具有不同的非平衡程度(IR)，IR等于負類數量正類數量.數據集名稱、數據量、特征維度數以及非平衡程度如表1所示.

在分類準確性評價上，本文采用AUC指標作為評價標準.AUC是公認的對非平衡數據分類評價較為合適的指標，AUC越大表示分類器的性能越好[14].為了評價算法的概率預測準確性，采用Brier指標和標度損失來評價分類結果，它們的定義如下：

Table 1 Description of Datasets Used in Our Experiments表1 真實數據集概況

3.2 形狀參數求取實驗

在形狀參數尋優實驗中，測試的方法是最大化后驗概率和貝葉斯分析.2種方法的初始ξ值均為0.1.圖2展示最大化后驗概率在9組生成數據集下ξ值隨迭代次數的變化圖像.

從圖2可以看到，9組生成數據集下的ξ值在最大化后驗概率方法下基本收斂，且ξ值變化的方向與其正確值一致，這說明最大化后驗概率方法對于形狀參數尋優的有效性.最終搜尋到的ξ值如表2所示.

(a) N=100, real ξ=-0.1 (b) N=500, real ξ=-0.1 (c) N=1000, real ξ=-0.1 (d) N=100, real ξ=0.2 (e) N=500, real ξ=0.2 (f) N=1000, real ξ=0.2 (g) N=100, real ξ=0.5 (h) N=500, real ξ=0.5 (i) N=1000, real ξ=0.5

Fig. 2 Convergence of ξ under MAP method圖2 最大化后驗概率方法下ξ值的收斂情況

Table 2 Estimation of the Shape Parameter Under MAP Method

觀察表2可知，最大化后驗概率方法最終搜尋到的ξ值基本正確，但當數據量為100時ξ值離正確值最遠，數據量為1 000時最近.這表明最大化后驗概率的參數尋優性能隨數據量的增大而提高，數據量過少時ξ值有較大偏差，但仍在可以接受的范圍內.

接下來是貝葉斯分析對參數尋優的實驗.實驗中，EM算法的迭代次數為10，每次迭代調用MH算法采樣的樣本數量為10 000～20 000均勻遞增.圖3為貝葉斯分析下對數似然概率值以及ξ值的收斂過程.

從圖3可以看到，采樣方法下的參數收斂過程并不平滑，細微的抖動情況非常多.ξ值基本朝正確的方向收斂，表示目標函數優化正確.貝葉斯分析方法最終搜尋到的ξ值如表3所示.

(a) Log Probability when N=100, real ξ=-0.1 (b) Log Probability when N=500, real ξ=-0.1 (c) Log Probability when N=1000, real ξ=-0.1

(d) ξ value when N=100, real ξ=-0.1(e) ξ value when N=500, real ξ=-0.1(f) ξ value when N=1000, real ξ=-0.1(g) Log Probability when N=100, real ξ=0.2(j) ξ value when N=100, real ξ=0.2(h) Log Probability when N=500, real ξ=0.2(k) ξ value when N=500, real ξ=0.2(i) Log Probability when N=1000, real ξ=0.2(l) ξ value when N=1000, real ξ=0.2(m) Log Probability when N=100, real ξ=0.5(p) ξ value when N=100, real ξ=0.5(n) Log Probability when N=500, real ξ=0.5(q) ξ value when N=500, real ξ=0.5(o) Log Probability when N=1000, real ξ=0.5(r) ξ value when N=1000, real ξ=0.5

Fig. 3 Convergence of ξ and Log Likelihood under Bayesian analysis圖3 貝葉斯分析方法下的ξ值和對數似然的收斂情況

Table 3 Estimation of the Shape Parameter Under Bayesian Analysis

從最終參數尋優結果來看，數據量為100時貝葉斯分析有一個徹底錯誤的數據，其他結果基本正確.還可以發現數據量的多少對結果的影響沒有像最大化后驗概率方法那樣明顯，數據量為500和1 000的參數尋優結果相差并不大.

3.3 對比實驗

本節實驗測試GEV線性回歸算法在真實數據集上的表現，包括分類準確性和概率預測準確性.

對比的其他算法為Logit回歸、Probit回歸和Cloglog回歸算法，還包括用SMOTE技術預處理數據的版本[16]，其中用SMOTE技術預處理后數據的非平衡度大約為1.數據集隨機劃分為70%的訓練集和30%的測試集，實驗結果取10次分類的均值.本文的GEV線性回歸算法使用最大化后驗概率方法來選取最優的形狀參數ξ，形狀參數尋優的迭代次數和線性回歸的迭代次數分別為20次和100次.各算法在各數據集上的AUC指標、Brier指標和標度損失分別如表4～6所示，其中最優結果用粗體表示.所使用的統計顯著性檢驗方法為雙側Wilcoxon檢驗，結果后的符號“*”和“**”分別表示在置信水平為90%和95%時該結果與該行中的最好結果相比在統計上是顯著的.

在分類準確性上，GEV線性回歸算法在大部分真實數據集上都能取得較好的分類結果，而由Logit和Probit鏈接函數這2種對稱函數構成的回歸算法在SMOTE預處理技術的幫助下也能得到部分最優結果，這從另一個角度說明它們的對稱性只適用于平衡數據.在概率預測準確性上，GEV線性回歸算法給出了大部分的最優結果，表明其在分類的概率預測值上能提供更準確的結果.其他3種回歸算法在SMOTE預處理下并沒有給出更好的結果，很多情況下反而更糟糕，這可能是因為SMOTE預處理下新生成的人為數據擾亂了概率值的準確預測.

Table 4 Results of Classification Accuracy on the Datasets in Terms of AUC Score表4 以AUC指標評價分類準確度的結果

Notes： The superscripts “*” and “**” indicate statistical significance of the best method in that row with respect to other methods on the 90% and 95% confidence interval respectively.

Table 5 Results of Class Probability Estimation on the Datasets in Terms of Brier Score表5 以Brier指標評價概率預測準確度的結果

Continued (Table 5)

Notes： The superscripts “*” and “**” indicate statistical significance of the best method in that row with respect to other methods on the 90% and 95% confidence interval respectively.

Table 6 Results of Class Probability Estimation on the Datasets in Terms of Calibration Loss表6 以標度損失評價概率預測準確度的結果

Notes： The superscripts “*” and “**” indicate statistical significance of the best method in that row with respect to other methods on the 90% and 95% confidence interval respectively.

3.4 運行效率實驗

本節實驗內容主要有2部分:1)觀察在3種最優化方法(利普希茨優化(Lipschitz)、梯度下降法(Gradient)、牛頓法(Hessian))下目標優化函數的收斂速度，也即其值隨迭代次數的增多而下降的趨勢;2)簡單測量3種最優化方法的運行時間.其中梯度下降法還混合采用了二分線搜索方法加快步長參數的搜索.

1) 對3個最優化方法的收斂速度的實驗.在作圖觀察時，由于校準損失函數無法求值，所以用對數損失函數代替.這里選取了6個具有代表性的真實數據集，各取了前10輪迭代，結果如圖4所示:

Fig. 4 Comparison of convergence rate of three optimization methods圖4 3種最優化方法的收斂速度比較

從圖4可以看出，隨著算法不斷迭代，3種最優化方法都能使對數損失函數變小，其中梯度下降法減小的速度最慢，利普希茨優化次之，牛頓法最快，這符合預期.另外還可以看到，利普希茨優化的收斂速度其實與牛頓法非常接近，兩曲線幾乎重合.這得益于利普希茨優化在保留了數據的二階結構的同時用常數近似了海森矩陣，取得了與牛頓法差不多的效率.

2) 方法運行時間的實驗.3種最優化方法在每個真實數據集上都運行了10遍，每種最優化方法也都跑滿了100次迭代.最終結果如表7所示，用時最短的結果已加粗顯示.顯然，利普希茨優化取得了所有真實數據集的最優結果，這得益于其將海森矩陣的逆矩陣計算提前到了迭代之前，大大減少了每輪迭代的運算量.可以看到，本文使用的利普希茨優化在算法運行效率上有顯著的提高.

Table 7 Comparison of Average Running Time of Three Optimization Methods

4 小 結

為了解決非平衡數據分類的概率預測問題，本文提出了GEV分布作為鏈接函數與校準損失函數相結合的算法.具有形狀參數的GEV分布可以提供高度靈活的曲線形狀來適應各種非平衡數據，并且目標優化函數是一個凸函數，算法確保了全局最優解的存在.利用GEV逆鏈接函數的利普希茨連續性，本文提出的GEV線性回歸算法在運行效率上得到了顯著改善.此外，本文還提出了2個基于概率的模型來求取形狀參數，并利用人工數據集上的實驗證實了2個模型對形狀參數求取的有效性.

在真實數據集上的實驗顯示了本文算法有著較高的概率預測準確性以及良好的分類性能.此外，在與其他優化方法的比較中，本文所提優化方法表現出很高的計算效率.未來的工作包括繼續研究更健壯的形狀參數求取方法以及對高維數據的實驗.