單純形分布聯(lián)合位置與散度模型的貝葉斯變量選擇

趙遠英，段星德，龐一成

（1.貴陽學(xué)院 數(shù)學(xué)與信息科學(xué)學(xué)院，貴陽 550005；2.楚雄師范學(xué)院 數(shù)學(xué)與統(tǒng)計學(xué)院，云南 楚雄 675000；3.貴州財經(jīng)大學(xué) 數(shù)學(xué)與統(tǒng)計學(xué)院，貴陽 550025）

0 引言

單純形(simplex)分布是一類非常重要的統(tǒng)計分布，Barndorff-Nielsen和Jorgensen[1]給出了其最早的分布形式。由于該分布是定義在區(qū)間(0,1)上的連續(xù)分布函數(shù)，在實際應(yīng)用時是分析百分比數(shù)據(jù)的有力統(tǒng)計工具，而且在公共衛(wèi)生、環(huán)境和經(jīng)濟等領(lǐng)域都有廣泛應(yīng)用，所以對單純形回歸模型的研究引起眾多統(tǒng)計工作者的普遍關(guān)注。例如：Song和Tan[2]基于廣義估計方程的方法研究了散度參數(shù)為常數(shù)的單純形分布廣義線性模型；之后Song等[3]基于Song和Tan[2]的方法對可變散度參數(shù)單純形分布廣義線性模型進行進一步的研究；Zhang和Wei[4]與Zhao等[5]得到單純形分布非線性混合效應(yīng)模型的極大似然估計與貝葉斯估計等。由于均值（位置）結(jié)構(gòu)與方差（散度）結(jié)構(gòu)建模的科學(xué)意義，其研究也成為統(tǒng)計研究的熱點之一。例如：Smyth[6]獲得了可變散度廣義線性模型的極大似然估計；基于聯(lián)合均值（位置）與方差（散度）結(jié)構(gòu)建模的思想，Paula[7]、Lee和Nelder[8]提出了雙重廣義線性模型的相關(guān)統(tǒng)計方法。但是上述文獻都未涉及到變量選擇問題。經(jīng)典的貝葉斯變量選擇方法有BIC準則[9]、貝葉斯因子[10]和偽貝葉斯因子[11]等。Mitchell和Beauchamp[12]在線性回歸模型的框架下提出一種貝葉斯子集選擇方法；George和McCulloch[13]基于隨機搜尋和Gibbs抽樣[14]有效識別潛在的解釋變量；Green[15]討論如何應(yīng)用可逆跳躍MCMC（馬爾科夫鏈蒙特卡羅）算法來研究模型的不確定性；Nott和Leonte基于Ising模型和Swendsen-Wang算法[16]描述了在廣義線性模型框架下的貝葉斯變量選擇抽樣算法。Leng等[17]利用貝葉斯自適應(yīng)LASSO方法研究如何選擇變量。盡管已經(jīng)有許多學(xué)者研究貝葉斯變量選擇方法，也做出了卓有成效的工作。然而，幾乎沒有文獻研究單純形分布聯(lián)合位置與散度模型的貝葉斯選擇問題。因此，本文的主要目的就是在貝葉斯統(tǒng)計和單純形分布聯(lián)合位置與散度模型的框架下，提出一種可行有效的變量選擇方法。

1 模型與符號

為了進行貝葉斯變量選擇，受Kuo和Mallick[18]思想的啟發(fā),本文考慮單純形分布聯(lián)合位置與散度模型：

其中，yi～S(μi，)表示響應(yīng)變量yi服從位置參數(shù)為μi∈(0，1)、散度參數(shù)為的單純形分布,其概率密度函數(shù)為：

2 貝葉斯變量選擇方法

記θ=(βT，αT)T，為了進行貝葉斯統(tǒng)計分析,首先得指定參數(shù)θ和示性變量向量γ與ω的先驗分布。假設(shè)θ，γ與ω的先驗分布相關(guān)獨立,且γj(j=1，2，…，p)，ωk(k=1，2，…，q)之間的先驗分布也相互獨立。即{θ，γ，ω}的先驗分布可以表示為:類似于Kuo和Mallick[18]與Zhao等[5]，本文考慮參數(shù)

β，α 和 γj(j=1，2，…，p) 與ωk(k=1，2，…，q)的如下先驗分布:

其中N(β0，Σβ)表示均值為β0協(xié)方差陣為 Σβ的多元正態(tài)分布,B(1，pj)表示參數(shù)為1和pj的Bernoulli分布，β0，Σβ，α0，Σα，和pj與ppk是事先給定的超參數(shù)。

本文提出一種能同時有效識別位置模型與散度模型中預(yù)測變量的貝葉斯變量選擇算法。記y=(y1，…，yn)T，x={x1，…，xn}和z={z1，…，zn} 。對參數(shù)θ，γ與ω的貝葉斯統(tǒng)計分析可以基于聯(lián)合后驗分布p(θ，γ，ω|y，x，z)∝p(y，x，z|θ，γ，ω)p(θ，γ，ω)， 其中p(y，x，z|θ，γ，ω) 是觀測數(shù)據(jù){y，x，z}的似然函數(shù)。因為分層模型的復(fù)雜性,很容易發(fā)現(xiàn)直接從聯(lián)合后驗分布p(θ，γ，ω|y，x，z)中抽樣是不可能的。本文采用近代統(tǒng)計計算工具MCMC方法來解決上述困難。即采用Gibbs抽樣[11]方法得到聯(lián)合后驗分布p(θ，γ，ω|y，x，z)的隨機樣本序列{(θ(t)，γ(t)，ω(t)):t=1，2，…，T}，然后基于此隨機序列,通過 γj(j=1，2，…，p)與ωk(k=1，2，…，q)的后驗概率來有效識別預(yù)測變量xij和zik，最終達到貝葉斯變量選擇的目的。從聯(lián)合后驗分布p(θ，γ，ω|y，x，z)中產(chǎn)生的隨機樣本序列的Gibbs抽樣步驟為：（1）選取參數(shù)θ，γ與ω的初始值為θ(0)，γ(0)和ω(0),并令t=0 ；（2）從條件分布p(θ||y，x，z，γ(t)，ω(t))中抽取θ(t+1)；（3）從條件分布p(γ||y，x，z，θ(t+1)，ω(t))中抽取γ(t+1)；（4）從條件分布p(ω||y，x，z，θ(t+1)，γ(t+1))中抽取ω(t+1)；（5）令t=t+1，重復(fù)步驟（2）至步驟（4）直到算法收斂。易知,θ，γ與ω的聯(lián)合后驗分布p(θ，γ，ω|y，x，z)正比于γ，ω)， 即：

2.1 條件分布

由于θ=(βT，αT)T,本文分別給出β和α的條件分布。

對β的條件分布，由式（3）、式（4）和式（5）有：

對α的條件分布，由式（3）、式（4）和式（5）有：

對γ和ω的條件分布,由式（3）、式（4）和式（5）分別有：

2.2 完成抽樣

以下簡述如何通過式（6）至式（9）產(chǎn)生所需要的隨機樣本。首先描述如何通過式（8）和式（9）產(chǎn)生所需要的γ和ω。對任何的j=1，2，…，p，記γ-j為γ刪除第j分量之后得到的向量。在給定y，x，z，β，α，γ-j，ω的條件下，可以通過Bernoulli分布向量γ*(γ**)的第j分量等于1(0)，其余分量和γ相等；相似地，對任何的k=1，2，…，q，記ω-k為ω刪除第k分量之后得到的向量。在給定y，x，z，β，α，γ，ω-k的條件下,也可以通過Bernoulli分布ω**)，向量ω*(ω**)的第k分量等于1(0)，其余分量和ω相等。

由上面的描述可知,從條件分布式（8）和式（9）中抽樣是很容易的。但是式（6）與式（7）中的條件分布并不是常見的標準分布,因此從條件分布式（6）與式（7）中抽樣就十分困難和復(fù)雜。本文采用Metropolis-Hastings(MH)算法[19,20]來解決這一困難。對于式（6）中條件分布的抽樣，本文選取多元正態(tài)分布N(0，Ωβ)為建議分布,其中

從條件分布中p(β|y，x，z，α，γ，ω)抽β的MH算法的步驟如下：已知β在第t次迭代時的迭代值為β(t),從分布N(β(t)，Ωβ)中隨機抽取潛在的轉(zhuǎn)移點 β*,并同時獨立地從區(qū)間[0,1]上的均勻分布U(0，1)中抽取一個隨機數(shù)u，如令β(t+1)=β(t)。類似地，對于式（7）中條件分布的抽樣，選取多元正態(tài)分布為建議分布，其中Ωα=(Vα+

從條件分布中p(α|y，x，z，β，γ，ω)抽α的MH算法的步驟如下：已知α在第t次迭代時的迭代值為α(t),從分布N(α(t)，Ωα)中隨機抽取潛在的轉(zhuǎn)移點 α*,并同時獨立地從區(qū)間[0,1]上的均勻分布U(0，1)中抽取一個隨機數(shù)τ，如否則令α(t+1)=α(t)。在具體應(yīng)用時,通常推薦選取σβ2與σα2使得接受β與α的平均概率為0.25左右。

3 數(shù)值例子

假設(shè)響應(yīng)變量yi(i=1，…，n)來自式（1）中定義的單純形分布聯(lián)合位置與散度模型,位置參數(shù)μi(i=1，…，n)與散度參數(shù)σi2(i=1，…，n)分別由下面的模型確定：

考慮p=4與q=3的變量選擇問題，假設(shè)預(yù)測變量xi=(xi1，xi2，xi3，xi4)T與zi=(zi1，zi2，zi3)T分別獨立地從標準正態(tài)分布N(0，I4)與N(0，I3)中產(chǎn)生，其中Im表示m×m的單位矩陣。參數(shù)的真值被設(shè)置為:(β1γ1，β2γ2，β3γ3，β4γ4)T=(0，0，1，1)T， (α1ω1，α2ω2，α3ω3)T=(0，0，1)T。即換言之，式（10）與式（11）分別變?yōu)?

因此，式（1）、式（10）和式（11）定義的模型為包含模型不確定性的單純形分布聯(lián)合位置與散度模型的類型。在模擬研究時,數(shù)據(jù)由式（1）、式（12）和式（13）產(chǎn)生,考慮樣本量分別為n=50，100和200三種情形，然后應(yīng)用本文介紹的Gibbs抽樣和MH算法進行貝葉斯變量選擇分析。表1基于20000次Gibbs抽樣給出了出現(xiàn)頻率最高的3個模型的比例。由模擬研究的數(shù)值結(jié)果可以得到以下結(jié)論：（1）貝葉斯變量選擇方法能同時正確地識別出位置模型中的預(yù)測變量為xi3和xi4，散度模型中的預(yù)測變量是zi3；（2）隨著樣本量n的增大，貝葉斯變量選擇方法的表現(xiàn)越來越好。

表1 模擬研究的數(shù)值結(jié)果

4 結(jié)論

本文針對單純形分布聯(lián)合位置與散度模型，基于Gibbs抽樣與MH算法進行貝葉斯變量選擇研究。結(jié)果表明，所提出來的方法能同時對位置模型與散度模型進行變量選擇。模擬研究的結(jié)果顯示了模型與方法的可行性與有效性。