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

趙遠英，段星德，龐一成

（1.貴陽學院 數學與信息科學學院，貴陽 550005；2.楚雄師范學院 數學與統計學院，云南 楚雄 675000；3.貴州財經大學 數學與統計學院，貴陽 550025）

0 引言

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

1 模型與符號

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

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

2 貝葉斯變量選擇方法

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

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

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

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

2.1 條件分布

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

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

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

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

2.2 完成抽樣

以下簡述如何通過式（6）至式（9）產生所需要的隨機樣本。首先描述如何通過式（8）和式（9）產生所需要的γ和ω。對任何的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）中抽樣就十分困難和復雜。本文采用Metropolis-Hastings(MH)算法[19,20]來解決這一困難。對于式（6）中條件分布的抽樣，本文選取多元正態分布N(0，Ωβ)為建議分布,其中

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

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

3 數值例子

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

考慮p=4與q=3的變量選擇問題，假設預測變量xi=(xi1，xi2，xi3，xi4)T與zi=(zi1，zi2，zi3)T分別獨立地從標準正態分布N(0，I4)與N(0，I3)中產生，其中Im表示m×m的單位矩陣。參數的真值被設置為:(β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）分別變為:

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

表1 模擬研究的數值結果

4 結論

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