海量數(shù)據(jù)下廣義線性模型參數(shù)的聚合估計(jì)算法研究

陳少東，李志強(qiáng)

(北京化工大學(xué) 數(shù)理學(xué)院，北京 100029)

一、引言

Nelder和Wedderburn等在研究指數(shù)族分布時(shí)引入廣義線性模型[1]。該模型涵蓋了許多常用的統(tǒng)計(jì)模型，如：經(jīng)典線性模型、對數(shù)線性模型、Logistic模型和Probit模型等。廣義線性模型既可用于刻畫連續(xù)型數(shù)據(jù)，也可用于刻畫計(jì)數(shù)數(shù)據(jù)、屬性數(shù)據(jù)等離散型數(shù)據(jù)，因此廣泛應(yīng)用于生物學(xué)、醫(yī)學(xué)、經(jīng)濟(jì)學(xué)以及社會科學(xué)等領(lǐng)域。

廣義線性模型的參數(shù)估計(jì)方法及其性質(zhì)已經(jīng)得到了深入地研究。然而，隨著數(shù)據(jù)采集技術(shù)的迅猛發(fā)展，可用于分析的數(shù)據(jù)呈現(xiàn)海量化趨勢。傳統(tǒng)的估計(jì)方法在解決海量數(shù)據(jù)下模型的參數(shù)估計(jì)問題時(shí)常常會遇到計(jì)算機(jī)存儲空間不足等問題，因此，有必要探究新的算法來解決海量數(shù)據(jù)下廣義線性模型的參數(shù)估計(jì)問題。

分治算法是當(dāng)前解決海量數(shù)據(jù)下統(tǒng)計(jì)模型參數(shù)估計(jì)問題的主流方法之一，該算法采用分而治之的思想，首先將原始海量數(shù)據(jù)集劃分為若干方便處理的小數(shù)據(jù)集，然后在每塊小數(shù)據(jù)集上構(gòu)造統(tǒng)計(jì)量將數(shù)據(jù)進(jìn)行壓縮，最后構(gòu)造一個(gè)聚合算法將每塊小數(shù)據(jù)集上的統(tǒng)計(jì)量聚合成海量數(shù)據(jù)集上的估計(jì)結(jié)果。分治算法可以解決模型實(shí)現(xiàn)過程中將所有數(shù)據(jù)一次性加載到內(nèi)存中而帶來的內(nèi)存溢出問題。由于在每塊小數(shù)據(jù)集的數(shù)據(jù)壓縮過程中，數(shù)據(jù)塊之間是獨(dú)立進(jìn)行的，所以該算法不僅適用于單機(jī)實(shí)現(xiàn)，也適用于并行實(shí)現(xiàn)。在實(shí)現(xiàn)時(shí)采用常用的海量數(shù)據(jù)處理工具Hadoop、Spark即可輕松實(shí)現(xiàn)并行。目前已有一些基于Hadoop和Spark的算法研究[2-4]，這些研究結(jié)果都表明集群計(jì)算能有效提高計(jì)算效率。

分治算法已廣泛應(yīng)用于許多海量數(shù)據(jù)統(tǒng)計(jì)模型研究。其中，方方等研究了海量數(shù)據(jù)下基于分治算法的線性模型和廣義線性模型的模型平均算法[5]；Lin和Xi等提出了一種基于分治算法的海量數(shù)據(jù)非線性估計(jì)方程的估計(jì)思想，其處理步驟為[6]：首先將海量數(shù)據(jù)集劃分為K塊子集，然后在每塊子集上基于估計(jì)方程將原始數(shù)據(jù)進(jìn)行壓縮，最后通過構(gòu)造聚合估計(jì)方程把每塊數(shù)據(jù)的壓縮信息有效的結(jié)合起來，從而得到模型系數(shù)的聚合估計(jì)。

然而，針對海量數(shù)據(jù)廣義線性模型的估計(jì)及其統(tǒng)計(jì)推斷還有待于進(jìn)一步深入研究。事實(shí)上，方方等[5]提出的平均估計(jì)難以進(jìn)行有效的統(tǒng)計(jì)推斷；而Lin和Xi僅研究了具有典則連接函數(shù)時(shí)聚合估計(jì)結(jié)果的漸近性質(zhì)，并不適用于一般的廣義線性模型。因此海量數(shù)據(jù)下具有一般連接函數(shù)的廣義線性模型的估計(jì)算法與統(tǒng)計(jì)推斷還有待于進(jìn)一步研究。

本文將Lin和Xi提出的非線性聚合估計(jì)方法推廣到具有一般連接函數(shù)的廣義線性模型的估計(jì)算法之中，提出了模型系數(shù)的聚合擬似然估計(jì)算法。首先，將整個(gè)數(shù)據(jù)集劃分為單機(jī)可處理的子數(shù)據(jù)集；其次，基于擬似然估計(jì)方程和Newton-Raphson算法的公式，將每塊子數(shù)據(jù)集壓縮為部分低階統(tǒng)計(jì)量進(jìn)行保存；然后根據(jù)保存的統(tǒng)計(jì)量和逼近的擬似然函數(shù)可構(gòu)造出聚合估計(jì)方程；最后求解得到聚合擬似然估計(jì)。該方法避免了在估計(jì)時(shí)使用全部數(shù)據(jù)進(jìn)行迭代計(jì)算，從而有效地克服了傳統(tǒng)估計(jì)方法在海量數(shù)據(jù)下帶來的計(jì)算機(jī)存儲空間不足等問題。由于聚合擬似然估計(jì)方程把每塊子數(shù)據(jù)的壓縮信息有效的結(jié)合起來，所以得到的聚合估計(jì)要比僅對每塊數(shù)據(jù)的估計(jì)結(jié)果做簡單的算術(shù)平均更能有效地利用樣本數(shù)據(jù)信息。而且根據(jù)定理1可以看出，當(dāng)子數(shù)據(jù)集的數(shù)量控制在一定范圍內(nèi)時(shí)，聚合擬似然估計(jì)與基于全部數(shù)據(jù)得到的擬似然估計(jì)是漸近等價(jià)的，因此可以得到良好的統(tǒng)計(jì)性質(zhì)。模擬結(jié)果表明，聚合擬似然估計(jì)算法不僅在解決存儲空間不足的同時(shí)提高了計(jì)算效率，而且有效地應(yīng)用于集群計(jì)算中。

二、廣義線性模型及其估計(jì)方法

若隨機(jī)變量Y的概率分布或概率密度函數(shù)滿足如下形式：

其中a(·)，b(·)，c(·)為已知函數(shù)，θ稱為典則參數(shù)(Canonical Parameter)，φ稱為擴(kuò)散參數(shù)(Dispersion Parameter)，則稱Y服從指數(shù)族分布。

設(shè)(Yi，Xi)，Xi∈Rp，i=1，2，…，N，為來自(Y，X)的獨(dú)立樣本，對于給定的X=xi，若滿足：

(1)條件密度函數(shù)f(yi|X=xi)服從指數(shù)族分布；

(2)設(shè)μi=E(yi|X=xi)為yi的數(shù)學(xué)期望，若存在已知的單調(diào)、可微函數(shù)g(·)滿足：

這里β∈Rp為未知參數(shù)，g稱為連接函數(shù)，則稱(Y，X)服從廣義線性模型[7]。

在廣義線性模型的建模過程中，常用的連接函數(shù)有[7]：

(1)恒等函數(shù)

g(μ)=μ，(μ∈R)

(2)對數(shù)函數(shù)

g(μ)=lnμ，(μ>0)

(3)Logistic函數(shù)

g(μ)=ln(μ/(1-μ))，(0<μ<1)

(4)Probit函數(shù)

g(μ)=Φ-1(μ)，(0<μ<1)

(5)重對數(shù)函數(shù)

g(μ)=ln(-lnμ)，(μ>0)

(6)互補(bǔ)重對數(shù)函數(shù)

g(μ)=ln(-ln(1-μ))，(0<μ<1)

其中，Φ為標(biāo)準(zhǔn)正態(tài)分布的累積函數(shù)。針對某個(gè)具體的廣義線性模型，絕大部分都是非典則連接函數(shù)，因此其在海量數(shù)據(jù)下的估計(jì)方法與統(tǒng)計(jì)推斷更具有一般性。

為了估計(jì)模型的未知參數(shù)β，將對數(shù)似然函數(shù)關(guān)于β求導(dǎo)可得到估計(jì)方程：

(1)

利用以上表達(dá)式，當(dāng)Y的指數(shù)族分布未知時(shí)，通常采用以下的擬似然估計(jì)方程方法。

設(shè)Q(μ，x)為擬似然函數(shù)[7]，滿足條件：

=GTV-1(Y-g*)=0

(2)

(3)

Newton-Raphson算法在求解過程中，每一步迭代需使用全部觀測數(shù)據(jù)，在觀測數(shù)N很大時(shí)，可能存在內(nèi)存不足等問題，因此有必要在海量數(shù)據(jù)情形下對該算法進(jìn)行改進(jìn)。

三、海量數(shù)據(jù)下的擬似然方程估計(jì)算法

(一)AQMLE算法

現(xiàn)給出聚合估計(jì)方程方法的思想。首先基于全部數(shù)據(jù)的擬似然估計(jì)方程為：

(4)

(5)

(6)

(7)

(8)

(二)AQMLE算法實(shí)現(xiàn)

AQMLE算法的計(jì)算量主要在數(shù)據(jù)壓縮步驟，而對于每一塊子數(shù)據(jù)集，因?yàn)樵跀?shù)據(jù)壓縮步驟沒有數(shù)據(jù)塊與數(shù)據(jù)塊之間的通信問題，所以很適合通過并行算法實(shí)現(xiàn)。在海量數(shù)據(jù)下，Hadoop與Spark是常用的處理工具，這些工具以集群的方式來完成海量數(shù)據(jù)的計(jì)算問題。

圖1 并行計(jì)算流程

四、AQMLE估計(jì)的漸近性質(zhì)

條件1.均值函數(shù)g-1(·)三階連續(xù)可導(dǎo)且一階導(dǎo)不為零；

條件2.方差函數(shù)V(·)>0且二階連續(xù)可導(dǎo)；

條件4.x為固定設(shè)計(jì)且一致有界；

(9)

(10)

注：當(dāng)廣義線性模型(1)中的連接函數(shù)為典則連接函數(shù)，根據(jù)指數(shù)族分布和廣義線性模型的定義可知，模型的條件數(shù)學(xué)期望、方差函數(shù)滿足以下條件：

(11)

通過直接計(jì)算可證：

因此擬似然估計(jì)方程(4)可簡化為：

(12)

結(jié)合式(11)，通過直接計(jì)算可得：

再由式(12)可知，對于具有典則連接函數(shù)的廣義線性模型，本文定理2的結(jié)論可簡化為與Lin和Xi定理5.2一致的結(jié)論。

五、數(shù)值模擬

本節(jié)對AQMLE算法進(jìn)行模擬檢驗(yàn)，以檢驗(yàn)算法的有效性。模擬分為兩個(gè)部分，第一部分是在普通單機(jī)環(huán)境下，AQMLE算法的結(jié)果分析；第二部分是在普通單機(jī)組成的Spark集群下，AQMLE算法的結(jié)果分析。有效性將從時(shí)間消耗與估計(jì)精度兩個(gè)角度來衡量。模擬時(shí)，我們選擇了常用的兩點(diǎn)分布，設(shè)yi|xi～b(1，pi)，其中pi可以通過以下兩種方式來建模：

(1)在典則連接下，有Logistic模型：

i=1，2，…，N

(2)在非典則連接下，選擇Probit模型：

其中Φ(·)為標(biāo)準(zhǔn)正態(tài)分布的累積函數(shù)。下面用本文所提出估計(jì)方法分別對上述兩個(gè)模型進(jìn)行模擬計(jì)算，并通過計(jì)算結(jié)果分析估計(jì)方法的有效性。

(一)單機(jī)計(jì)算

在單機(jī)情況下，取數(shù)據(jù)量N=3×106，β0的維數(shù)為32，分塊數(shù)K在1到300之間，自變量獨(dú)立同分布，服從標(biāo)準(zhǔn)正態(tài)分布。本次實(shí)驗(yàn)在Ubuntu16.0.4操作系統(tǒng)下使用Python3.7.0完成，所有實(shí)驗(yàn)重復(fù)50次取平均值作為最終結(jié)果。

表1(a)和表1(b)給出的是算法在不同分塊數(shù)K下的計(jì)算時(shí)間和估計(jì)偏差結(jié)果。首先關(guān)注計(jì)算效率部分，從表中可以看出，兩個(gè)模型計(jì)算消耗的時(shí)間都隨著K的增大先下降后趨于穩(wěn)定。當(dāng)K=1時(shí)，計(jì)算消耗的時(shí)間最多，在K=10時(shí)，計(jì)算速度有顯著提升，而當(dāng)K繼續(xù)增大時(shí)，由于分塊計(jì)算本身會帶來程序的開銷，所以計(jì)算效率不再上升，會逐漸趨于穩(wěn)定狀態(tài)。在Logistic模型中，算法將計(jì)算效率提升了30%左右，而Probit模型的計(jì)算效率提升了20%左右。

表1(a) Logistic模型在單機(jī)環(huán)境下取不同K的試驗(yàn)結(jié)果

表1(b) Probit模型在單機(jī)環(huán)境下取不同K的試驗(yàn)結(jié)果

表2 不同K下隨機(jī)抽樣估計(jì)優(yōu)勢比例表

(二)集群計(jì)算

在實(shí)際中，面對GB甚至是TB級別的數(shù)據(jù)量時(shí)，數(shù)據(jù)可能是以分布式的方式存儲在多臺計(jì)算機(jī)上的。幸運(yùn)的是，本文的AQMLE方法在分布式框架上的實(shí)現(xiàn)是簡單的。這部分我們將在Spark計(jì)算框架下進(jìn)行模擬。

本次實(shí)驗(yàn)使用的模型及數(shù)據(jù)規(guī)模仍與單機(jī)計(jì)算時(shí)一致。實(shí)驗(yàn)在Ubuntu16.0.4下通過4臺普通計(jì)算機(jī)組成的Spark集群來完成，計(jì)算環(huán)境為Python3.7.0、Hadoop2.7.7及Spark2.3.2。其中，計(jì)算節(jié)點(diǎn)的配置為3.40GHzCPU、4G內(nèi)存、8核。

表3展示的是Spark集群下的計(jì)算結(jié)果，該表給出了不同分塊數(shù)下Logistic模型和Probit模型的計(jì)算時(shí)間。如表所示，兩個(gè)模型的實(shí)驗(yàn)結(jié)果大致相同：計(jì)算消耗的時(shí)間先有一個(gè)下降趨勢，然后趨于穩(wěn)定。當(dāng)分塊數(shù)K取值過小時(shí)，由于部分計(jì)算資源被閑置，會導(dǎo)致計(jì)算時(shí)間略高。但是，之后隨著計(jì)算資源的充分利用，計(jì)算時(shí)間逐漸下降，此時(shí)集群的計(jì)算效率明顯優(yōu)于單機(jī)環(huán)境。

表3 集群環(huán)境下實(shí)驗(yàn)結(jié)果

六、實(shí)證分析

本節(jié)，我們使用來自UCI數(shù)據(jù)庫的SUSY數(shù)據(jù)集(http：//archive.ics.uci.edu/ml/datasets/SUSY)建立Probit分類器，以此來檢驗(yàn)AQMLE算法在實(shí)際應(yīng)用中的可行性和有效性。SUSY數(shù)據(jù)集是一個(gè)超對稱粒子數(shù)據(jù)集，在高能物理領(lǐng)域非常受歡迎，由P.Baldi等首次使用[9]。該數(shù)據(jù)集共包含500萬個(gè)樣本，每個(gè)樣本有18個(gè)特征及1個(gè)標(biāo)簽，具體數(shù)據(jù)可通過UCI數(shù)據(jù)庫下載。該分類任務(wù)是為了判斷通過撞機(jī)碰撞后產(chǎn)生的粒子是否是感興趣的粒子。數(shù)據(jù)中前8個(gè)特征由加速器中的粒子測量儀測量得到，后10個(gè)特征是物理學(xué)家構(gòu)建的用于提高分類準(zhǔn)確率的組合特征。

本次實(shí)證分析分別基于前8個(gè)特征(記為lo-level)、后10個(gè)特征(記為hi-level)以及所有特征(記為lo+hi-level)構(gòu)建Probit分類器，并用前450萬條數(shù)據(jù)做訓(xùn)練集，后50萬條數(shù)據(jù)做測試集。我們選擇機(jī)器學(xué)習(xí)常用的ROC(Receiver Operating Characteristic)曲線來刻畫得到的Probit分類器的性能，并用AUC值(Area Under the Curve)作為評價(jià)指標(biāo)來進(jìn)行結(jié)果分析。通過對K取K=1，20，…，100，200進(jìn)行比較，發(fā)現(xiàn)基于上述三組特征構(gòu)建的Probit分類器的分類效果并不受K取值的影響，同一組特征構(gòu)建的Probit分類器在不同K值下的AUC值是一樣的，結(jié)果如表4示。該結(jié)果進(jìn)一步表明AQMLE算法在處理海量數(shù)據(jù)的問題上是可行。

表4 不同特征下分類器的AUC值

圖2 ROC曲線圖

從圖2中我們可以看到僅使用lo-level特征訓(xùn)練的分類器比僅使用hi-level特征訓(xùn)練的分類器的效果要差得多。然而，僅使用hi-level特征訓(xùn)練的分類器比使用lo+hi-level特征訓(xùn)練的分類器具有更弱的性能。這表明了基于lo-level特征構(gòu)建的Probit分類器沒有足夠的能力識別標(biāo)簽，而基于hi-level特征構(gòu)建的Probit分類器很大程度上提高了識別標(biāo)簽的能力，但hi-level特征也丟失了lo-level特征中部分可用于解釋標(biāo)簽的信息，該實(shí)驗(yàn)符合Baldi P等(2014)中給出的結(jié)論。

七、結(jié)論

本文提出了一種基于分治算法和Newton-Raphson迭代算法的方法用于處理海量數(shù)據(jù)下的擬似然估計(jì)問題，目的在于解決海量數(shù)據(jù)參數(shù)估計(jì)帶來的存儲問題和計(jì)算效率問題。首先，將整個(gè)數(shù)據(jù)集隨機(jī)劃分成K塊，基于Newton-Raphson算法將每塊子數(shù)據(jù)集通過構(gòu)造統(tǒng)計(jì)量的方式進(jìn)行壓縮并保存；然后基于保存的統(tǒng)計(jì)量導(dǎo)出整個(gè)數(shù)據(jù)集的聚合估計(jì)，記為AQMLE。分治算法符合并行計(jì)算的思想，因此AQMLE算法適用于超大規(guī)模的并行化問題。本文證明了當(dāng)K以一定速度趨于無窮時(shí)，基于AQMLE算法得到的估計(jì)結(jié)果仍具有漸近正態(tài)性。模擬與實(shí)證分析表明，AQMLE算法在解決計(jì)算過程中的存儲問題的同時(shí)提高了計(jì)算效率，且估計(jì)結(jié)果要優(yōu)于隨機(jī)抽樣和平均方法，而基于集群的實(shí)現(xiàn)方法可進(jìn)一步提高計(jì)算效率，因此，本文的工作可為海量數(shù)據(jù)下的模型參數(shù)估計(jì)問題提供一些參考。