基于偏正態(tài)數(shù)據(jù)下位置、均值回歸模型的參數(shù)估計

王丹璐，吳劉倉，鄭桂芬

(昆明理工大學(xué)理學(xué)院，云南昆明650504)

1.引言

生活中常見的統(tǒng)計分布可分為兩類:對稱分布和偏態(tài)分布.我們所學(xué)的大部分分布如正態(tài)分布、t分布、Laplace分布等都屬于對稱分布，故對稱分布在統(tǒng)計理論與方法的研究中占有重要地位.但是現(xiàn)實生活中收集到的大部分?jǐn)?shù)據(jù)都不具有嚴(yán)格的對稱性，幾乎都存在偏斜，這時再用對稱分布描述數(shù)據(jù)的性質(zhì)并不科學(xué)，所以引入了偏態(tài)分布來解決這一問題.偏態(tài)分布具有偏斜、尖峰、厚尾的特征，由于偏態(tài)分布比對稱分布更能捕捉到全面準(zhǔn)確、及時有效的信息，所以近年來偏態(tài)分布引起了許多學(xué)者的研究興趣，如鄭等[1]研究了基于偏Laplace正態(tài)數(shù)據(jù)下的位置、均值回歸模型的參數(shù)估計，朱等[2]研究了偏t正態(tài)數(shù)據(jù)下混合線性聯(lián)合位置與尺度模型的參數(shù)估計.此外，由于偏態(tài)分布在數(shù)據(jù)擬合運(yùn)用方面更廣闊，所以在農(nóng)業(yè)、環(huán)境、生物醫(yī)學(xué)應(yīng)用及金融經(jīng)濟(jì)等領(lǐng)域也有很多研究意義.

偏正態(tài)分布作為正態(tài)分布的推廣，不僅具有正態(tài)分布的良好統(tǒng)計特性，同時也有具有偏態(tài)分布的特征，適用性更廣泛，所以國內(nèi)外許多學(xué)者研究了偏正態(tài)分布問題.自Azzalini[3]首次定義了偏正態(tài)分布的概念，隨后關(guān)于它的性質(zhì)、參數(shù)估計、統(tǒng)計診斷、推廣以及多元偏正態(tài)分布已經(jīng)被系統(tǒng)深入地研究.本文中我們使用的是Sahu等[4]提出的多元偏正態(tài)分布的單變量形式，該偏態(tài)分布是基于偏橢圓分布變換而來的，相比Azzalini提出的偏正態(tài)分布該函數(shù)更具有吸引力，因為該分布對于偏度參數(shù)的刻畫較簡單，使用EM算法更容易實現(xiàn).在Sahu提出的偏正態(tài)分布模型中已有部分研究，比如Cancho等[5]在經(jīng)典的貝葉斯方法的基礎(chǔ)上，對偏正態(tài)非線性回歸模型進(jìn)行了深入的探討;Bolfarine等[6]提出了多元偏正態(tài)線性模型的四種診斷措施;Lachos等[7]在多元偏正態(tài)回歸模型基于極大似然估計的EM算法，研究了多元偏正態(tài)回歸模型的參數(shù)估計;Arellano-Valle等[8]研究了偏正態(tài)混合線性模型;Arellano-Valle等[9]研究了其貝葉斯推斷.

目前的對于該偏正態(tài)分布中，并沒有對位置和均值建模的內(nèi)容，但考慮到位置和均值在偏態(tài)分布中的重要地位，所以本文詳細(xì)介紹了利用EM算法和梯度下降法對位置回歸模型和均值回歸模型的參數(shù)估計，并通過模擬和實例結(jié)果表明本文提出的模型是可行的.

本文的目的是討論偏正態(tài)位置、均值回歸模型的參數(shù)估計，全文結(jié)構(gòu)安排如下:第二部分介紹了偏正態(tài)分布的部分性質(zhì);第三部分分別介紹了偏正態(tài)分布下位置和均值回歸模型;第四部分利用梯度下降法和EM算法分別對位置和均值回歸模型的參數(shù)進(jìn)行了估計;第五部分通過Monte Carlo模擬證實本文提出方法的有效性;第六部分通過實例分析并且和多元線性回歸模型比較，表明本文研究的模型和方法是科學(xué)合理的.

2.偏正態(tài)分布

根據(jù)Sahu(2003)等提出偏正態(tài)函數(shù)，若隨機(jī)變量Y服從偏正態(tài)分布，那么它的概率密度函數(shù)(pdf)記做:

其中φ(·|a，b2)和Φ(·|a，b2)分別表示N(a，b2)的密度函數(shù)(pdf)和分布函數(shù)(cdf)，記做Y～SN(μ，σ2，λ)，其中μ為位置參數(shù)，σ2為尺度參數(shù)，λ為偏度參數(shù).當(dāng)λ＜0時，曲線較長的尾部在左邊，稱作負(fù)偏分布;當(dāng)λ＞0時，曲線較長的尾部在右邊，稱作正偏分布;當(dāng)λ=0時，偏正態(tài)分布重新退化為對稱的正態(tài)分布.

I偏正態(tài)分布下的隨機(jī)表示

設(shè)X0～N(0，1)，X1～N(μ，σ2)，是兩個獨(dú)立隨機(jī)變量，則隨機(jī)變量Y～SN(μ，σ2，λ)的隨機(jī)表示為

從而得到偏正態(tài)分布的層次表示

即該偏正態(tài)分布可以分層為一個條件正態(tài)分布Y|(T=t)和一個半正態(tài)分布因此(Y，T)的聯(lián)合密度函數(shù)為

根據(jù)貝葉斯準(zhǔn)則有

3.位置和均值回歸模型

I位置回歸模型

由偏正態(tài)分布的概率密度函數(shù)(2.1)和位置回歸參數(shù)模型(3.1)可以得到

其中yi為第i個響應(yīng)變量，位置參數(shù)為μi，尺度參數(shù)為σ2，偏度參數(shù)為λ，xi=(xi1，xi2，···，xip)T為解釋變量，β=(β1，β2，···，βp)T為p×1維的位置回歸模型的未知參數(shù).

II均值回歸模型

其中α=(α1，α2，···，αp)T為p×1維的均值回歸模型的未知參數(shù).

因為有潛變量的存在，所以利用普通極大似然估計方法對參數(shù)估計比較困難，故我們使用基于梯度下降法的EM算法來解決含有潛變量問題的參數(shù)估計.接下來研究模型參數(shù)的極大似然估計的EM算法.

4.極大似然估計的EM算法和梯度下降法

自Dempster等[10]提出EM算法(Expectation Maximization Algorithm)來處理不完全數(shù)據(jù)以來，EM算法及其擴(kuò)展應(yīng)用已經(jīng)被廣泛的研究和應(yīng)用到各個領(lǐng)域.EM算法是一種常用的迭代算法，每一次迭代由兩部分組成:E-step，M-step.E-step是根據(jù)參數(shù)初始值或上一次迭代所得到的結(jié)果來計算出對數(shù)似然函數(shù)的期望值;M-step是通過E-step計算得到的對數(shù)似然函數(shù)的期望值，求出該似然函數(shù)的最大的參數(shù)值，記為新的參數(shù)值，用該參數(shù)值代替初始值或上一次迭代所得到的結(jié)果使得對數(shù)似然函數(shù)最大化.E-step和M-step兩步不斷交替計算，直到滿足收斂條件時停止迭代，從而得到的參數(shù)估計值更加的逼近真實的參數(shù)值.

下面給出EM算法在偏正態(tài)數(shù)據(jù)下位置回歸模型和均值回歸模型的參數(shù)估計中的具體計算步驟:

I位置回歸模型下極大似然估計的EM算法

由(3.2)式可得偏正態(tài)數(shù)據(jù)下位置回歸模型的似然函數(shù)為

令θ1=(βT，σ2，λ)T，則L(βT，σ2，λ)=L(θ1)，兩邊取自然對數(shù)，得到對數(shù)似然函數(shù)

設(shè)T為潛變量，利用位置回歸模型和偏正態(tài)分布的層次表達(dá)法得

其中c1為獨(dú)立于θ1的常數(shù).同時，因為該似然函數(shù)是單峰可導(dǎo)的，則似然方程的解存在唯一，且為極大似然估計，因而是強(qiáng)相合的.為了得到θ1的極大似然估計，將(4.3)式極大化，但此時極大化得到的估計依賴于潛變量.因此，為了處理潛變量，引入EM算法，使用觀測數(shù)據(jù)yi求出完全數(shù)據(jù)下對數(shù)似然函數(shù)的條件期望.

E-step:利用k次迭代，得到估計

利用(2.5)和(2.6)式給出的條件期望，可以計算出(4.4)式中的條件期望E和有

M-step:兩步長梯度法在構(gòu)造搜索方向時，充分利用上次迭代點信息，使計算不易陷入局部最大值的陷阱，故結(jié)合梯度下降法[11]將對數(shù)似然函數(shù)極大化，獲取新的參數(shù)值，給出兩個初值設(shè)計以下迭代

其中

利用(4.7)式給出的目標(biāo)函數(shù)，得到位置參數(shù)的Score函數(shù)為

利用梯度下降法將E-step和M-step反復(fù)迭代，直至算法收斂.

II均值回歸模型下極大似然估計的EM算法

令θ2=(αT，σ2，λ)T，則L(αT，σ2，λ)=L(θ2)，兩邊取自然對數(shù)，得到對數(shù)似然函數(shù)

利用均值回歸模型和偏正態(tài)分布的層次表達(dá)法得

用層次表達(dá)法得到完全數(shù)據(jù)下的對數(shù)似然函數(shù)為

其中c2為獨(dú)立于θ2的常數(shù).

E-step:利用k次迭代，得到估計:

利用(2.5)和(2.6)式給出的條件期望，可以計算出(4.10)的條件期望和有

皮膚表面脂質(zhì)是由甘油三酯、蠟酯、角鯊烯、脂肪酸以及少量膽固醇、膽固醇酯和雙甘酯組成的非極性脂類混合物[14]，其中角鯊烯具有高度的不飽和性，在UV誘導(dǎo)下，角鯊烯脂質(zhì)過氧化形成角鯊烯過氧化物（SQOOH）以及丙二醛（MDA）[15]，增加羰基化反應(yīng)、糖基化反應(yīng)，從而導(dǎo)致皮膚泛黃。

M-step:利用梯度下降法將對數(shù)似然函數(shù)極大化，獲取新的參數(shù)值，給出兩個初值設(shè)計以下迭代

其中

利用(4.14)式給出的目標(biāo)函數(shù)，得到均值參數(shù)的Score函數(shù)為

利用梯度下降法將E-step和M-step反復(fù)迭代，直至算法收斂.

III極大似然估計的漸近性質(zhì)

同時考慮該極大似然估計的漸近正態(tài)性，取θ0為θ的真值，假定以下正則條件

(i)xi=(xi1，xi2，···，xip)T是固定的;

(ii)參數(shù)空間是緊的，真值θ0為參數(shù)空間的內(nèi)點;

(iii)yi，i=1，2，···，n相互獨(dú)立，Y～SN(μ，σ2，λ)，其中在位置回歸下的極大似然估計時，;在均值回歸下的極大似然估計時，.

在假定(i)-(iii)之下，似然方程L(θ)=0在L(1)(θ)=0在n→+∞時有相合解滿足且似然方程任一相合解必為θn的BAN估計.

5.Monte Carlo模擬

I位置回歸模型下的Monte Carlo模擬

為評價上述位置回歸模型參數(shù)估計方法的有效性，使用該位置回歸模型參數(shù)估計方法對有限樣本進(jìn)行Monte Carlo模擬，并用均方誤差(MSE)、偏差(Bias)來評價該參數(shù)估計的精確度，其定義如下

其中θ1=(βT，σ2，λ)，θ1(0)是參數(shù)θ1的真值.根據(jù)模型(3.1)，考慮偏正態(tài)數(shù)據(jù)下位置參數(shù)的回歸模型.

根據(jù)模型(5.1)產(chǎn)生模擬數(shù)據(jù)，其中解釋變量x的分量獨(dú)立產(chǎn)生于均勻分布U(-1，1)，其余參數(shù)的真值設(shè)定為-2.

II均值回歸模型下的Monte Carlo模擬

為評價上述均值回歸模型參數(shù)估計方法的有效性，使用該位置均值模型參數(shù)估計方法對有限樣本進(jìn)行Monte Carlo模擬，并用均方誤差(MSE)、偏差(Bias)來評價該參數(shù)估計的精確度，其定義如下

其中θ2=(αT，σ2，λ)，θ2(0)是參數(shù)θ2的真值.根據(jù)模型(3.3)，考慮偏正態(tài)數(shù)據(jù)下均值參數(shù)的回歸模型.

分別取樣本量n=100，200，300，400，重復(fù)模擬1000次.模擬結(jié)果見表1，評價結(jié)果見表2、表3.

表1 位置、均值回歸模型的參數(shù)估計模擬結(jié)果

表2 位置回歸模型的參數(shù)估計評價結(jié)果

表3 均值回歸模型的參數(shù)估計評價結(jié)果

從上表可以得到，無論樣本是左偏還是右偏，即參數(shù)λ是負(fù)值還是正值，隨著樣本量n的增大，所有參數(shù)的估計值越來越接近真值，而且通過估計量的均方誤差(MSE)和偏差(Bias)的判斷也可以得到當(dāng)樣本量不斷增大，估計效果越好.以上結(jié)論表明，本文提出的偏正態(tài)數(shù)據(jù)下位置和均值回歸模型及所使用的EM算法對參數(shù)的極大似然估計取得了較理想的效果.

6.實例分析

根據(jù)2017年國際糖尿病聯(lián)名協(xié)會公布的一組數(shù)據(jù)，目前全球有4.7億糖尿病患者，其中中國糖尿病患者高達(dá)1.14億，位居世界第一，所以對糖尿病患者數(shù)據(jù)的研究是必要的.本文選取UCI機(jī)器學(xué)習(xí)數(shù)據(jù)庫的一組糖尿病數(shù)據(jù)進(jìn)行分析，數(shù)據(jù)包括了478名糖尿病患者的身體質(zhì)量指數(shù)BMI和其他相關(guān)指標(biāo)，該數(shù)據(jù)中包含一個響應(yīng)變量Y-糖尿病患者BMI，7個解釋變量X，分別是懷孕次數(shù)、血糖(mmol/L)、血壓(mmHg)、皮質(zhì)厚度(mm)、胰島素(pmol/ml)、糖尿病遺傳函數(shù)、年齡.利用R軟件計算得y偏度為0.5529，峰度為3.0376，從峰度值和偏度值，可以知道，糖尿病患者的BMI數(shù)據(jù)具有明顯的偏斜，而正態(tài)分布的偏度值為0，峰度值為3，直方圖如下:

圖1 糖尿病患者BMI的直方圖

根據(jù)圖1和偏度系數(shù)說明糖尿病患者的身體質(zhì)量指數(shù)BMI數(shù)據(jù)近似的服從偏正態(tài)分布，所以可以利用該組數(shù)據(jù)做偏正態(tài)分布的位置和均值回歸模型的參數(shù)估計，考慮Y與X之間的模型如下:

同時，我們加入多元線性回歸模型做比較，多元線性回歸模型如下:

糖尿病患者BMI的位置、均值回歸模型參數(shù)估計和多元線性回歸參數(shù)估計結(jié)果:

表4 糖尿病患者BMI的位置、均值回歸模型參數(shù)估計和多元線性回歸參數(shù)估計結(jié)果

由于在同一組糖尿病患者BMI的數(shù)據(jù)中，尺度參數(shù)和偏度參數(shù)應(yīng)為一致，從表格中分析到，位置回歸模型和均值回歸模型的尺度參數(shù)σ2和偏度參數(shù)λ基本一致，但是多元線性回歸模型的尺度參數(shù)σ2與其他兩個模型差別較大，并且不存在偏度參數(shù)λ.這是由于使用多元線性回歸模型需要忽略數(shù)據(jù)的偏斜特征.從解釋變量的系數(shù)來看，多元線性回歸模型中的懷孕次數(shù)γ1、血糖γ2、血壓γ3、皮脂厚度γ4、年齡γ7都不顯著，并且系數(shù)值與其余兩個模型相差過大，胰島素γ5呈現(xiàn)正相關(guān)，這都與實際情況不符.以上，我們可以判斷出該數(shù)據(jù)用多元線性回歸模型做參數(shù)估計導(dǎo)致估計效果不合理甚至錯誤.這是由于這組數(shù)據(jù)中并不是嚴(yán)格對稱的，具有一定的偏斜，但多元線性回歸模型針對的是對稱分布數(shù)據(jù)進(jìn)行統(tǒng)計推.所以在數(shù)據(jù)呈現(xiàn)偏正態(tài)情形下，不應(yīng)該再選用多元線性回歸模型做參數(shù)估計，可以選擇本文提出的偏正態(tài)數(shù)據(jù)下的位置回歸模型和均值回歸模型做數(shù)據(jù)分析.

在均值和位置回歸模型中，解釋變量中血糖X2、皮脂厚度X4、糖尿病遺傳函數(shù)X6呈現(xiàn)正相關(guān)且系數(shù)較大，這與實際情況相符合.血糖升高是判斷是否患糖尿病重要指標(biāo)，肥胖是二型糖尿病的主要誘因，并且如果父母患糖尿病那么子女患糖尿病的風(fēng)險比正常人要高.解釋變量中胰島素X5呈負(fù)相關(guān)，這與患有糖尿病的患者缺少胰島素實際情況相符合.

7.結(jié)論

Monte Carlo模擬表明了本文提出的EM算法對位置和均值回歸模型的未知參數(shù)進(jìn)行了較好的估計，與現(xiàn)有的模型和常用估計方法相比較，本文提出的偏正態(tài)數(shù)據(jù)下位置和均值回歸模型具有更好的運(yùn)用性，更適合實際生活中的偏態(tài)數(shù)據(jù)的研究.除此以外，在實例分析中，對糖尿病患者BMI數(shù)據(jù)的應(yīng)用研究也表明了本文提出的模型和方法能更準(zhǔn)確分析解釋數(shù)據(jù)，因此表明本文提出的模型和方法是可行的.