朱志娥,吳劉倉(cāng),戴 琳

(昆明理工大學(xué)理學(xué)院,云南昆明650093)

偏t正態(tài)數(shù)據(jù)下混合線性聯(lián)合位置與尺度模型的參數(shù)估計(jì)

朱志娥,吳劉倉(cāng),戴 琳

(昆明理工大學(xué)理學(xué)院,云南昆明650093)

偏t正態(tài)分布是分析尖峰,厚尾數(shù)據(jù)的重要統(tǒng)計(jì)工具之一.研究提出了偏t正態(tài)數(shù)據(jù)下混合線性聯(lián)合位置與尺度模型,通過EM算法和Newton-Raphson方法研究了該模型參數(shù)的極大似然估計(jì).并通過隨機(jī)模擬試驗(yàn)驗(yàn)證了所提出方法的有效性.最后,結(jié)合實(shí)際數(shù)據(jù)驗(yàn)證了該模型和方法具有實(shí)用性和可行性.

偏t正態(tài)分布;混合線性聯(lián)合位置與尺度模型;EM算法;極大似然估計(jì)

§1 引 言

在社會(huì)經(jīng)濟(jì)領(lǐng)域中,經(jīng)常收集到的數(shù)據(jù)存在著大量分類問題且不是嚴(yán)格地服從正態(tài)分布,而是服從具有明顯偏斜的厚尾分布.面對(duì)這些情況,如果簡(jiǎn)單的對(duì)數(shù)據(jù)總體進(jìn)行分析,那么就很難找到數(shù)據(jù)中各類別的差異.因此,為了更好的分析數(shù)據(jù)的特征,常常對(duì)數(shù)據(jù)進(jìn)行聚類分析.聚類分析是將總體數(shù)據(jù)按不同指標(biāo)或特性進(jìn)行分類,再對(duì)具有相似性質(zhì)或相似指標(biāo)的數(shù)據(jù)進(jìn)行詳細(xì)的分析,這使得在研究大數(shù)據(jù)時(shí)減少了很多工作量.一方面,偏t正態(tài)分布能夠較好的刻畫具有明顯偏斜和厚尾數(shù)據(jù)的分布規(guī)律.另一方面,了解是哪些因素影響了數(shù)據(jù)出現(xiàn)偏斜和厚尾的情況,這一問題也是研究者們所關(guān)注的熱點(diǎn).因此,對(duì)偏t正態(tài)數(shù)據(jù)的統(tǒng)計(jì)分析具有重要的理論價(jià)值和實(shí)際意義,本文主要研究偏t正態(tài)數(shù)據(jù)下混合線性聯(lián)合位置與尺度模型的參數(shù)估計(jì)問題.

混合回歸模型是研究含有兩個(gè)及兩個(gè)以上子聚類的混合數(shù)據(jù)的最重要的統(tǒng)計(jì)分析工具之一,在經(jīng)濟(jì)學(xué)、醫(yī)學(xué)、環(huán)境科學(xué)、生物學(xué)、抽樣調(diào)查及工程技術(shù)等領(lǐng)域具有廣泛的應(yīng)用.自Goldfeld和Quandt[1]首次介紹了混合回歸模型以來,混合回歸模型就吸引了很多研究學(xué)者的關(guān)注,可參見McLachlan和Peel[2]的專著及其里面的參考文獻(xiàn).最近,Yao等[3]基于t分布研究提出了穩(wěn)健的混合線性回歸模型;Song等[4]基于Laplace分布研究提出了穩(wěn)健的混合線性回歸模型;Azzalini[5]系統(tǒng)地對(duì)偏正態(tài)分布的性質(zhì)進(jìn)行了研究;Kotz和Vicari[6]綜述了偏態(tài)分布的研究現(xiàn)狀;Nadarajah等[7]也討論了一系列偏態(tài)分布及其各自的性質(zhì),并求得了這些分布的特征函數(shù)和n階矩的表達(dá);Liu和Lin[8]基于偏正態(tài)分布,首次研究了偏態(tài)數(shù)據(jù)下混合線性回歸模型.

以上所研究的數(shù)據(jù)特征僅為偏斜或厚尾,并未同時(shí)對(duì)具有偏斜和厚尾的數(shù)據(jù)進(jìn)行研究,且以上研究都假定混合數(shù)據(jù)的每個(gè)子聚類的方差都相同,偏度也相同.然而,在許多實(shí)際問題中這樣的假定是不合理的.針對(duì)異方差數(shù)據(jù),很多統(tǒng)計(jì)學(xué)者對(duì)聯(lián)合均值與方差模型進(jìn)行了系統(tǒng)的研究.如吳劉倉(cāng)等[9]研究了偏t正態(tài)數(shù)據(jù)下聯(lián)合位置與尺度模型的極大似然估計(jì);Aitkin[10]給出了聯(lián)合均值與方差模型的極大似然估計(jì);吳劉倉(cāng)等[11]研究了聯(lián)合均值與方差模型的變量選擇.此外,G′omez[12]針對(duì)偏t正態(tài)模型討論了其極大似然估計(jì)方法和Fisher信息矩陣及其相關(guān)性質(zhì);吳劉倉(cāng)等[13]研究了缺失偏t正態(tài)數(shù)據(jù)下線性回歸模型的參數(shù)估計(jì);Lin等[14]研究了偏t正態(tài)數(shù)據(jù)下非線性回歸模型中方差的齊次性檢驗(yàn).然而,針對(duì)混合偏態(tài)異方差數(shù)據(jù)研究較少.

故本文針對(duì)偏t正態(tài)混合數(shù)據(jù)、異方差和線性回歸,提出了偏t正態(tài)數(shù)據(jù)下混合線性聯(lián)合位置與尺度模型.其次利用EM算法和Newton-Raphson方法研究此模型參數(shù)的極大似然估計(jì).最后,結(jié)合實(shí)際數(shù)據(jù)說明該模型和方法具有實(shí)用性和可行性.

本文由以下幾個(gè)部分構(gòu)成：§2首先簡(jiǎn)要的介紹偏t正態(tài)分布及其密度函數(shù)的表示,在其基礎(chǔ)上提出了偏t正態(tài)數(shù)據(jù)下混合線性聯(lián)合位置與尺度模型;§3詳細(xì)的介紹了偏t正態(tài)數(shù)據(jù)下混合線性聯(lián)合位置與尺度模型的EM算法,并給出了此模型下EM算法的操作步驟;§4通過隨機(jī)模擬試驗(yàn)驗(yàn)證了該方法的有效性;§5結(jié)合實(shí)際數(shù)據(jù)說明了該模型和方法的實(shí)用性和可行性;§6是本文的小結(jié)和討論.

§2 偏t正態(tài)數(shù)據(jù)下混合線性聯(lián)合位置與尺度模型

2.1 偏t正態(tài)分布

偏t正態(tài)分布實(shí)際是一種廣義的正態(tài)分布.考慮隨機(jī)變量Y服從StN分布,同G′omez[12]一樣,即可以表示成Y ～StN(μ,σ2,λ,ν),其中,μ為位置參數(shù),σ 為尺度參數(shù),λ為偏度參數(shù),ν為自由度,隨機(jī)變量Y的概率密度函數(shù)可表示為：

其中,tν(·)為自由度為ν的t分布的概率密度函數(shù),Φ(·)為標(biāo)準(zhǔn)正態(tài)分布的分布函數(shù).顯然,當(dāng)自由度ν→∞時(shí),概率密度函數(shù)(1)即為SN分布的概率密度函數(shù),當(dāng)偏度參數(shù)λ=0時(shí),將得到自由度為ν的t分布的概率密度函數(shù).當(dāng)自由度ν→∞和λ=0時(shí),該分布為正態(tài)分布.在本文中,只考慮自由度ν是有限的情況.當(dāng)自由度未知時(shí),可通過截面似然的方法來對(duì)其進(jìn)行估計(jì).在偏t正態(tài)分布中,有

2.2 偏t正態(tài)數(shù)據(jù)下混合線性聯(lián)合位置與尺度模型

設(shè)y1,y2,...,yn是n個(gè)來自m個(gè)混合比例為π1,π2,...,πm的隨機(jī)樣本量,概率密度函數(shù)為f(y)=StN(μ,σ2,λ,ν),即Y 服從混合偏t正態(tài)分布,

本文中,假定混合子聚類數(shù)m是固定和已知的.當(dāng)然,在一些實(shí)際應(yīng)用中,m可能是未知的,它也需要和混合比例以及其它參數(shù)一起估計(jì),但在本文中,為了簡(jiǎn)化，僅考慮是m已知的情形,只對(duì)未知參數(shù)θ進(jìn)行估計(jì).

本文針對(duì)偏t正態(tài)混合數(shù)據(jù),異方差和線性回歸,提出如下偏t正態(tài)數(shù)據(jù)下混合線性聯(lián)合位置與尺度模型：

在模型(3)中,{yi}1是獨(dú)立的響應(yīng)變量,xi={xi1,...,xip}T和hi={hi1,...,hiq}T是解釋變量,βj={βj1,...,βjp}T是第j個(gè)子聚類中維數(shù)為p×1的位置模型的未知參數(shù),γj={γj1,...,γjq}T是第j個(gè)子聚類中維數(shù)為q×1的尺度模型的未知參數(shù),λj是第j個(gè)子聚類的偏度參數(shù).兩個(gè)解釋變量可能完全相同,部分相同或者完全不同,即位置模型和尺度模型可能包含一些相同的解釋變量或者包含不同的解釋變量.即使包含相同的解釋變量也可能分別對(duì)位置模型和尺度模型產(chǎn)生不同的影響.

§3 偏t態(tài)數(shù)據(jù)下混合線性回歸模型的EM算法

Mclachlan和Peel[2]指出EM算法可以獲得有限混合模型任意分布的極大似然估計(jì).本文假定子聚類m是固定和已知的,由于每一次觀測(cè)值都屬于m個(gè)子聚類中的一類,但無(wú)法確定是m類中的第幾類.為了解決這個(gè)問題,引入潛變量zi=(zi1,...,zim),

對(duì)不可確定來自m類中的第幾類的情況,可以利用多項(xiàng)式分布zij進(jìn)行刻畫,若yi屬于m個(gè)子聚類中的第j類,則zij等于1,若yi不屬于m個(gè)子聚類中的第j類,而屬于m-1類中的某一類,此時(shí)zij等于0. 則完全數(shù)據(jù)下關(guān)于參數(shù)的對(duì)數(shù)似然函數(shù)可以寫成：

EM算法(Expectation Maximization Algorithm)是一種聚類算法,其具體流程分為兩個(gè)步驟進(jìn)行：E-step和M-step.第一步,E-step是根據(jù)參數(shù)初始值或上一次迭代所得結(jié)果來計(jì)算對(duì)數(shù)似然函數(shù)的期望值;第二步,M-step是將對(duì)數(shù)似然函數(shù)最大化以獲得新的參數(shù)值,用新得到的參數(shù)值代替初始值或上一次迭代所得結(jié)果使得對(duì)數(shù)似然函數(shù)最大化.重復(fù)執(zhí)行以上兩步驟,直至收斂.

下面給出EM算法在偏t正態(tài)數(shù)據(jù)下混合線性聯(lián)合位置與尺度模型中的計(jì)算步驟：E-step：利用從k次迭代得到的估計(jì)θ(k),計(jì)算出：

EM算法的操作步驟：

步驟2：E-step,在k+1次迭代中,通過引入潛變量zij來確定來自m類中那一類的數(shù)據(jù),然后計(jì)算完全數(shù)據(jù)下對(duì)數(shù)似然函數(shù)的條件期望：

步驟3：M-step,在k+1次迭代中,使用下面的公式進(jìn)行最大化的計(jì)算：

利用Newton-Raphson方法能得出θ的極大似然估計(jì).

步驟4：重復(fù)步驟2和步驟3,直至收斂.

當(dāng)自由度ν未知時(shí),可以通過最大化對(duì)數(shù)似然函數(shù)將自由度ν和混合回歸模型參數(shù)θ一起估計(jì)出來.為了簡(jiǎn)化,也可以考慮ν的截面最大似然(Pro fi le likelihood)估計(jì)：

紅松是國(guó)家保護(hù)植物，其物種的價(jià)值珍貴，工業(yè)，農(nóng)業(yè)，城市建筑和室內(nèi)裝飾。因此，植樹造林工程栽培受到重點(diǎn)關(guān)注。松樹造林工程保護(hù)自然資源和生態(tài)環(huán)境，但也產(chǎn)生經(jīng)濟(jì)效益的，造林工程項(xiàng)目，對(duì)造林工程技術(shù)的高要求。

注：一方面,因?yàn)槠玹正態(tài)分布是單峰的;另一方面,當(dāng)|θ(k)-θ(k-1)|≤δ,δ為預(yù)定的充分小的正數(shù),如δ=10-3等,保證了EM算法的收斂.

§4 Monte Carlo模擬

為了評(píng)價(jià)上述參數(shù)估計(jì)方法的有效性,對(duì)有限樣本性質(zhì)進(jìn)行模擬研究,參數(shù)估計(jì)的精確度使用均方誤差(MSE)來評(píng)價(jià)和衡量,其定義如下：

其中βj(0),γj(0),λj(0),πj(0)分別是βj,γj,λj,πj的真值.

根據(jù)下面的模型(10),產(chǎn)生模擬數(shù)據(jù)：

其中yi(i=1,2,...,n)是根據(jù)StN分布混合產(chǎn)生的相互獨(dú)立的響應(yīng)變量,且yi服從偏t正態(tài)分布.xi～U(-1,1),hi～U(-1,1)是各自分量相互獨(dú)立的解釋變量,分別取π1(0)=0.25和0.50由π2=1-π1可知：π2(0)=0.75和0.50給定β1,γ1,λ1和β2,γ2,λ2的真值分別為：β1(0)=(0,1,1)T,γ1(0)=(0,1,1)T,λ1(0)=0.5,β2(0)=(0,-1,-1)T,γ2(0)=(0,-1,-1)T,λ2(0)=-0.5 取樣本量n=400,600,800,分別取ν=3,4,5,重復(fù)模擬1000次.通過1000次模擬的均方誤差的平均數(shù)評(píng)價(jià)和衡量參數(shù)估計(jì)的精度,模擬結(jié)果見表1,表2和表3.根據(jù)表1,表2和表3的模擬結(jié)果,可得到,無(wú)論樣本是左偏或者右偏,即偏度參數(shù)λ取負(fù)值或者正值：

表1 ν=3的模擬結(jié)果

表2 ν=4的模擬結(jié)果

表3 ν=5的模擬結(jié)果

(1)在相同的自由度ν下,隨著樣本量n的增加所有參數(shù)的估計(jì)值越來越接近真值,且估計(jì)的均方誤差MSE也越來越小.

以上結(jié)論表明,本文基于偏t正態(tài)數(shù)據(jù)提出的混合線性回歸模型及所使用的EM算法對(duì)參數(shù)的極大似然估計(jì)取得了較理想的效果.

§5 實(shí)例分析

下面利用本文提出的偏t正態(tài)數(shù)據(jù)下混合線性聯(lián)合位置與尺度模型及其方法,對(duì)運(yùn)動(dòng)員體質(zhì)指數(shù)數(shù)據(jù)進(jìn)行統(tǒng)計(jì)推斷.

人體體質(zhì)指數(shù)(BMI)是衡量一個(gè)人是否肥胖或超重的標(biāo)準(zhǔn)之一.采用202名運(yùn)動(dòng)員的體質(zhì)指數(shù)(BMI)數(shù)據(jù)[16]進(jìn)行統(tǒng)計(jì)分析,其中女性運(yùn)動(dòng)員100名,男性運(yùn)動(dòng)員102名.體質(zhì)指數(shù)(BMI)包含有響應(yīng)變量Y-kg/(cm)2和8個(gè)解釋變量：X1-紅細(xì)胞計(jì)數(shù),X2-白細(xì)胞計(jì)數(shù),X3-比容,X4-血紅蛋白X5-等離子鐵蛋白濃度,X6-總皮膚褶皺,X7-體脂百分比,X8-去脂體重.將建立體質(zhì)指數(shù)Y(BMI)與8個(gè)解釋變量之間的關(guān)系.(這里假定xi與hi相同).

首先,對(duì)體質(zhì)指數(shù)(BMI)數(shù)據(jù)進(jìn)行正態(tài)性檢驗(yàn),并作出直方圖與QQ圖,見圖1,圖2及圖3,且分別為全部202名運(yùn)動(dòng)員,100名女性運(yùn)動(dòng)員和102名男性運(yùn)動(dòng)員的體質(zhì)指數(shù)(BMI)的正態(tài)檢驗(yàn)的概率圖.從以下這三幅圖可直觀上看出,數(shù)據(jù)近似服從正態(tài)分布,但具有明顯的偏斜,且均值附近與尾區(qū)的概率值比正態(tài)分布大,而其余區(qū)域的概率比正態(tài)分布小,即所研究數(shù)據(jù)具有尖峰、厚尾的特征.

圖1 202名運(yùn)動(dòng)員的體質(zhì)指數(shù)(BMI)的直方圖與QQ圖

圖2 100名女性運(yùn)動(dòng)員的體質(zhì)指數(shù)(BMI)的直方圖與QQ圖

利用模型(10),建立如下模型：

表4 體質(zhì)指數(shù)BMI參數(shù)的估計(jì)結(jié)果

從表4中可以得出：對(duì)于三組數(shù)據(jù)所建立的聯(lián)合位置模型和尺度模型有細(xì)微的差異.解釋變量X1-紅色細(xì)胞計(jì)數(shù),在聯(lián)合位置模型和尺度模型中對(duì)女性體質(zhì)指數(shù)的影響較小,對(duì)混合數(shù)據(jù)和男性數(shù)據(jù)的影響成一致關(guān)系;在兩個(gè)模型中解釋變量X3-比容,對(duì)混合數(shù)據(jù)和女性數(shù)據(jù)的影響成一致關(guān)系;在尺度模型中解釋變量X7-體脂百分比,對(duì)混合數(shù)據(jù)成正影響,而對(duì)男性數(shù)據(jù)和女性數(shù)據(jù)成負(fù)影響.若不考慮分類研究,其結(jié)果可能會(huì)偏離實(shí)際影響,從而得到錯(cuò)誤的結(jié)論.表4所得到的精細(xì)結(jié)果能為進(jìn)一步深入研究體質(zhì)指數(shù)(BMI)提供了重要的統(tǒng)計(jì)保障.

§6 結(jié)論

本文基于混合偏t正態(tài)數(shù)據(jù),異方差和線性回歸,研究提出了偏t正態(tài)數(shù)據(jù)下混合線性聯(lián)合位置與尺度模型,并利用EM算法和Newton-Raphson方法研究了該模型參數(shù)的極大似然估計(jì).Monte Carlo隨機(jī)模擬結(jié)果表明模型參數(shù)估計(jì)方法有效和可行.通過對(duì)運(yùn)動(dòng)員體質(zhì)指數(shù)(BMI)實(shí)際數(shù)據(jù)的分析也說明了該模型和方法具有實(shí)用性.

參考文獻(xiàn)：

[1] Goldfeld S M,Quandt R E.A Markov model for switching regression[J].Journal of Econometrics,1973,1：13-16.

[2] McLachlan G,Peel D.Finite Mixture Models[M].New York：Wiley，2000.

[3] Yao Wenxin,Wei Yan,Yu Chun.Robust mixture regression models using t-distribution[J].Computational Statistics and Data Analysis,2014,71：116-127.

[4] Song Weixing,Yao Wenxin,Xing Yanru.Robust mixture regression models fi tting by Laplace distribution[J].Computational Statistics and Data Analysis,2014,71：128-137.

[5] Azzalini A.A class of distributions which includes the normal ones[J].Scandinavian Journal of Statistics,1985,2：171-178.

[6] Kotz S,Vicari D.Survey of developments in the theory of continuous skewed distributions[J].International Journal of Statistics,2005,2：225-261.

[7] Nadarajah S,Kotz S.Skew distributions generated from di ff erent families[J].Acta Applicandae Mathematica,2006,91：1-37.

[8] Liu Min,Lin Tsung-I.A skew-normal mixture regression model[J].Educational and Psychological Measurement,2014,74：139-162.

[9] 吳劉倉(cāng),馬婷,戴琳.基于StN分布下聯(lián)合位置與尺度模型的極大似然估計(jì)[J].應(yīng)用數(shù)學(xué),2013,26：671-676.

[10]Aitkin M.Modelling variance heterogeneity in normal regression using GLIM[J].Applied Statistics,1987,36：332-339.

[11]吳劉倉(cāng),張忠占,徐登可.聯(lián)合均值與方差模型的變量選擇[J].系統(tǒng)工程理論與實(shí)踐,2012,32：1754-1760.

[12]G′omez H W,Venegas O,Bolfarine H.Skew-symmetric distributions generated by the distribution function of the normal distribution[J].Environmetric,2007,18：395-407.

[13] 吳劉倉(cāng),張家茂,李玲雪.缺失偏t正態(tài)數(shù)據(jù)下線性回歸模型的統(tǒng)計(jì)推斷[J].應(yīng)用數(shù)學(xué),2015,28：16-25.

[14]Lin Jinguan,Xie Fengchang,Wei Bocheng.Statistical diagnostics for skew-t-normal nonlinear models[J].Communications in Statistics-Simulation and Computation,2009,38：2096-2110.

[15]Taylor and Francis.Finite mixture of generalized semiparametric models variable slection via penalized estimation[J].Communications in Statistics-Simulation and Computation,2015,DOI：10.1080/03610918.2014.953687.

[16]Cook A D,Weisberg S.An introduction to regression graphics[M].John Wiley Sons,1994.

Parameter estimation for linear joint location and scale models with mixture skew-t-normal data

ZHU Zhi-e,WU Liu-cang,DAI Lin
(Faculty of Science,Kunming University of Science and Technology,Kunming 650093,China)

Skew-t-normal distribution is one of the most important statistical tools to analyze the obvious peak and fat tail data.A linear mixture joint location and scale model with skew-t-normal data is proposed in this paper.The maximum likelihood estimation of the unknown parameters of this model is investigated based on Expectation Maximization(EM)algorithm and Newton-Raphson method.Furthermore,the proposed procedure works satisfactorily through Monte Carlo experiments.Finally,a real example shows that both this model and method are useful and e ff ective.

skew-t-normal distribution；mixture of linear joint location and scale models；EM algorithm；maximum likelihood estimation

62F10;62J12

O212.1

A

：1000-4424(2016)04-0379-11

2016-03-24

2016-04-26

國(guó)家自然科學(xué)基金(11261025;11026309);云南省自然科學(xué)基金(2011FZ044)