一種二元響應變量模型的分布式貝葉斯估計方法

吳 磊, 錢夕元

(華東理工大學理學院,上海 200237)

一種二元響應變量模型的分布式貝葉斯估計方法

吳 磊, 錢夕元

(華東理工大學理學院,上海 200237)

在海量數據背景下,傳統的基于單個計算節點的算法很難滿足分析要求。考察了一種分布式貝葉斯估計方法,通過在每臺機器上單獨運行蒙特卡洛抽樣并做加權平均可以有效地解決算法效率問題。將該方法應用于基于廣義極值模型的二元響應變量回歸分析,并探討其實用性。模擬研究表明分布式算法比傳統方法更有效。

海量數據; 分布式貝葉斯方法; 極值模型

邏輯回歸模型是在處理二元響應變量數據時最為常用的一種廣義線性模型,它采用邏輯分布作為連接函數,可以實現利用連續型解釋變量來說明0-1二元響應變量的變化。該模型一般假設潛在變量的概率響應曲線是對稱的,即0-1二元響應變量中的各類樣本數目基本均衡,但當樣本數存在明顯不平衡時,邏輯回歸模型會嚴重違背對稱性的假設,帶來連接函數設定錯誤,使得模型參數估計存在較大的偏差和均方誤差[1]。

不平衡數據在實際應用中并不少見,它一般來源于某類稀有事件或現象發生概率較小的環境,且相關屬性的數據會具有明顯的偏度特征。為此,學者們提出了大量改進的連接函數用來靈活處理此類數據。最近,Wang等[2]提出了以廣義極值(GEV)分布作為連接函數的二元響應變量回歸模型,該連接函數比傳統GEV分布增加了一個形狀參數,新增的形狀參數不但沒有取值約束,而且可以更大幅度地調節偏度,使得該模型對非對稱或對稱的響應曲線都可以進行擬合,具有了更廣泛的靈活性,可以更好地處理二元不平衡數據。

值得注意的是模型的靈活性帶來了經典參數估計方法(極大似然估計)求解的困難性。隨著馬爾科夫鏈蒙特卡羅(MCMC)方法的發展,貝葉斯估計方法得到了更加廣泛的應用。貝葉斯估計方法可以有效地利用先驗信息,對小到中型樣本問題可以有效改善估計精度,但如何有效地將貝葉斯方法應用到海量數據分析成為近期一個研究熱點。隨著數據爆炸式增長,單個的處理器已經很難滿足人們的需求,一個可以想到的解決辦法是將數據分發到多個處理器上,但隨之帶來的問題就是如何解決各個節點間的信息交換,如何協調好各個處理過程,否則很容易出現死鎖或者串行化等問題。貝葉斯方法中廣泛采用的MCMC方法是基于馬氏鏈構造的,其當前狀態轉移的概率依賴于前一個狀態,這實際上和分布式的思想是有沖突的,因為馬氏鏈要求串行化而分布式要求的是并行化[3]。

本文首先給出了基于極值理論的二元響應變量回歸模型及其貝葉斯估計,針對海量數據應用環境,給出了分布式貝葉斯估計算法,設計模擬數據驗證了算法的有效性。

1 二元響應變量模型及其貝葉斯估計(Be)

Ψ(·)是一個取值在0和1之間的非負函數,標準情況下,設定Ψ(·)為累積分布函數,而稱Ψ-1(·)為連接函數。

常用二元響應變量模型的連接函數有logit連接Ψ-1(pi)=lg{pi/(1-pi)},probit連接Ψ-1(pi)=Φ-1(pi),Φ-1為標準正態分布的反函數,cloglog連接Ψ-1(pi)=-lg{-lg(pi)}等。上述連接函數都是對稱的,在處理不平衡數據時會出現較大的偏差和均方誤差。

Wang等[2]提出了基于廣義極值分布的二元響應變量回歸模型,其連接函數采用如下累積分布函數,具體表示為:

ξ為形狀參數,用以改變模型分布的偏度和尾部厚度。根據ξ的不同,該模型既可以表現出對稱性,也可以表現出非對稱性,可以很好地用來處理非平衡樣本數據[4]。本文采用基于MCMC算法的貝葉斯方法進行模型的參數估計,其Metropolis-Hastings(M-H)抽樣算法描述如下:

Step0選取待估參數的初始值θ0;

Step1從產生候選參數的密度函數q(θ(s-1),θ*)中獲得候選參數θ*;

Step2計算候選參數被接收的概率α(θ(s-1),θ*);

Step3以α(θ(s-1),θ*)的概率設θ(s)=θ*,或者以1-α(θ(s-1),θ*)的概率設θ(s)=θ(s-1);

Step4重復Step1,Step2和Step3S次;

Step5以θ0,θ1,…,θS作為參數θ的近似后驗密度,作均值估計得參數θ的點估計。

其中接受概率

式中,p(θ=θ*|y)表示后驗密度在點θ=θ*的取值,q(θ*;θ=θ(s-1))表示隨機變量θ的密度函數在θ=θ(s-1)處的取值。

2 分布式貝葉斯估計(DBe)方法

2.1模型參數估計

在海量數據背景下,上述估計方法將變得十分困難。本文提出借助分布式貝葉斯方法對模型參數進行估計。該方法的主要思想是根據現有的計算資源,在確保每份數據集的樣本容量足夠的情況下合理地把樣本數據拆分,為每份數據分配一個獨立的計算節點做蒙特卡洛抽樣,從而得到參數的貝葉斯后驗分布,最后根據一定的方式把每份數據的后驗分布整合成一個全局的后驗分布,其主要過程描述如下:

記y為全部的樣本數據(樣本數為n),ys是第s份數據(樣本數為ns),記θ為待估參數。假設數據集間相互獨立,則根據貝葉斯公式:

(1)

從式(1)中可以看到,每個部分的先驗分布變成了總體先驗的S次方根,這是為了保持整個系統中的先驗信息保持不變。另外,由于對先驗信息并不是很了解,本文采用了方差較大的無信息正態先驗。根據模型可知,似然函數為:

由于后驗分布比較復雜,本文采用自適應的正態隨機游走來逼近該目標后驗,其優勢是可以提高Metropolis-Hastings迭代過程的接受率,根據上一次的參數是否被接受來調整隨機游走的步長,從而自適應地保證更高的接受率[5]。

假定在第s個計算節點上得到了總共G個Metropolis-Hastings抽樣,以下述加權平均的方法得到總共S個計算節點的全局后驗分布參數估計:

式中,Ws稱為權重矩陣,一般可以取樣本方差的逆或在參數θ維數較高時作簡單的平均即可。DBe估計的步驟如下:

Step1 把樣本數據y分片為y1,…,ys;

Step2 由分離先驗信息p(θ)1/S,重復S次獨立的M-H抽樣(該部分算法步驟如上)得θsg～p(θ|ys),g=1,…,G;

Step3 對θsg,s=1,…,S;g=1,…,G加權,得到全局的后驗分布參數θg,g=1,…,G;

Step4 根據接收率α,過濾掉馬爾科夫鏈的前(1-α) G個參數;

2.2其他非參數估計策略

3 模擬研究

取樣本容量n=100 000,解釋變量取3個,均由服從均值為0、方差為1的正態隨機數生成,模型參數的真值為β0=1,β1=2,β2=3,β3=4,分布參數的真值ξ=2,產生模擬數據集50份。根據模擬數據集的樣本容量以及計算資源,設定S=12,G=1 250,α=0.8。

數據分布式處理環境如下:Hadoop2.x:1個Master節點和3個Slave節點,集群中CPU的數量Intel(R) Xeon(R) CPU E5-2620 0 @ 2.00 GHz*8*(1[M]+3[S])=32core,內存大小32 (1[M]+3[S])=128 GB,JVM參數設定-Xmx768 M,其中Container的大小為1 GB[6]。

表1給出了分布式貝葉斯方法和傳統分析方法的對比結果。從中可以看到,分布式貝葉斯方法優勢明顯,在本文的計算節點上,計算效率提升了近3倍。但由于計算時間取決于最慢的節點,且Hadoop環境下的Map階段(將數據拆分成12份)并沒有并行化而是僅使用了一個節點,因此計算效率并不呈線性提升。但是隨著數據量的進一步增大,分布式貝葉斯方法的優勢會更趨明顯。另外,在準確性和穩定性上,當M-H抽樣的初值為0,自適應隨機游走的步長分別是0.005和0.1的情況下,50次的重復實驗中分布式貝葉斯方法和傳統分析方法的參數估計精度差別不大,參見圖1的比較結果。

表1 算法對比分析結果Table 1 Comparative analysis of distributed Bayes vs.classical method

圖1 分布式貝葉斯方法和傳統方法參數 估計精度對比分析結果Fig.1 Comparative analysis of distributed Bayes vs.classical method

4 結束語

本文針對不平衡二元響應變量的海量數據,基于廣義極值分布進行回歸分析,借助于分布式貝葉斯方法(Metropolish-Hastings抽樣算法)進行參數估計,模擬研究表明該方法具有一定的計算優勢,應用于模擬數據分析中取得了較好的效果。今后可以進一步推廣到基于Gibbs抽樣的分布式貝葉斯分析中。

[1] CZADO C,SANTNER T J.The effect of link misspecification on binary regression inference [J].Journal of Statistical Planning and Inference,1992,33(2):213-231.

[2] WANG X,DEY D K.Generalized extreme value regression for binary response data:An application to B2B electronic payments system adoption [J].The Annals of Applied Statistics,2010,4(4):2000-2023.

[3] STEVEN S L,BLOCKER A W,BONASSI F V .Bayes and big data:The consensus Monte Carlo algorithm [J].International Journal of Management Science and Engineering Management,2016,11(2):78-88.

[4] GHOSH S K,MUKHOPADHYAY P,LU J C.Bayesian analysis of zero-inflated regression models [J].Journal of Statistical Planning and Inference,2006,136(4):1360-1375.

[5] CHIPMA H A,GEORGE E I,MCCULLOCH R E.BART:Bayesian additive regression trees [J].The Annals of Applied Statistics,2010,4(1):266-298.

[6] DEAN J,GHEMAWAT S.Mapreduce:Simplified data processing on large clusters [J].Communications of the ACM,2008,51(1):107-113.

ADistributedBayesianRegressionMethodforBinaryResponseMassiveData

WULei,QIANXi-yuan

(SchoolofScience,EastChinaUniversityofScienceandTechnology,Shanghai200237,China)

In the background of massive data,it is difficult to meet the analysis requirements for traditional one-node based algorithm.This paper considers a distributed Bayesian estimation method to solve the GEV based general linear regression model by running a separate Monte Carlo algorithm on each machine.The method is applied to regression analysis of binary response variables based on generalized extreme value model.The results show that the proposed distributed Bayesian regression algorithm is much faster than the traditional algorithm in the simulated data sets studying.

massive data; distributed Bayesian regression; GEV model

1006-3080(2017)04-0559-04

10.14135/j.cnki.1006-3080.2017.04.016

2016-10-31

國家高科技研究發展(“863”)計劃(2015AA20107);上海市經信委“軟件和集成電路產業發展專項資金”(140304)

吳 磊(1992-),男,上海人,碩士生,主要研究方向為統計計算。

錢夕元,E-mail:xyqian@ecust.edu.cn

TP301.6

A