非標準分布貝葉斯分析的WinBUGS實現*

徐州醫學院流行病與衛生統計教研室(221002) 曾 平 王 婷 何 鵬

WinBUGS是一套程序化貝葉斯統計分析軟件，應用者只需要指定數據的似然函數和參數的先驗分布，然后就可通過Gibbs抽樣得到參數的后驗樣本，這對沒有顯性表達式的貝葉斯后驗分布尤其重要。Win-BUGS軟件中提供了23種常見的分布用于指定似然函數和先驗分布〔1〕，能夠滿足大部分的貝葉斯統計分析。然而如果貝葉斯分析中的分布形式超出了Win-BUGS軟件所提供的范圍，那么需要應用一定的編程技巧和數學轉換以適于WinBUGS的程序要求。本文主要介紹非標準分布貝葉斯分析的WinBUGS實現技巧，這里所謂的非標準分布是相對于WinBUGS中的23種分布而言。

原理和方法

1．一個正態數據的WinBUGS程序

假設數據X=xi，i=1，…，n來源于獨立的正態分布f(xi|μ，τ)，目的在于通過X對參數μ，τ做出推斷。參數的似然函數指定獨立的無信息先驗 μ ～N(0，10－6)和 τ～gamma(0.001，0.001)，根據貝葉斯定理，μ，τ的后驗分布為

f(μ，τ|X)∝L(X|μ，τ)N(0，10－6)gamma(0.001，0.001)需要說明的是:(1)WinBUGS軟件中正態分布的參數為均數和精度(precision)〔1〕，即 τ =1/σ2;(2)μ，τ 的標準無信息獨立先驗應該是μ∝1和τ∝1/τ，都屬于不規則先驗(improper prior)〔2〕，但 WinBUGS 中不允許指定不規則先驗，因此分別指定為 N(0，10－6)和gamma(10－3，10－3)，此時各自的方差為 106和 103，如此大的方差反應了對參數信息的模糊認識，這種先驗又被稱作局部均勻先驗(locally uniform prior)或低信息先驗(low informative prior)〔3〕，避免出現不規則的后驗分布。對應的WinBUGS程序如下:

model{程序以 model開始所有程序都包含在“{}”中

for i in 1:n{

x［i］～dnorm(mu，tau)}3 ～4 行通過 for循環語句建立似然函數

mu～dnorm(0，1．0E-6)指定 μ 的先驗

tau ～dgamma(0．001，0．001)指定 τ的先驗

sigma＜-sqrt(1/tau)計算標準差

}

所有的參數或者通過計算得到，如標準差sigma由精度計算得到，或者必須指定一個分布，如對x和mu指定正態分布，對 tau指定 gamma分布，“＜-”的作用和“=”一樣，“～”表示“服從某個分布”。Win-BUGS中提供的23種常用分布，見表1。

表1 WinBUGS中的常用分布

2．非標準分布的WinBUGS程序上式可看作二項分布數據yi=1，i=1，…，n的似然函數，參數pi=exp(li－c)。常數c的存在等價于右邊每一項乘以因子exp(c)，顯然不會影響貝葉斯后驗推斷，c的作用在于保證pi在［0，1］之間，可選一個大的整數，如10000。WinBUGS程序如下:

c＜-10000設定常數c，在不合適時需要調整以保證 p∈［0，1］

for i in 1:n{

l［i］＜ -＊＊＊指定對數似然

p［i］＜ -exp(l［i］-c)計算二項分布參數 p

y［i］＜-1對所有y設定為1

y［i］～ dbern(p［i］)}設定二項分布似然(2)利用Poisson分布指定非標準分布

常數c不會影響貝葉斯后驗推斷。上式可看作參數r=1的負二項分布數據yi=0，i=1，…，n的似然函數，參數 pi=exp(li－c)，c的作用在于保證 pi在［0，1］之間，可選一個大的整數，如10000。WinBUGS程序如下:

c＜-10000設定常數c，在不合適時需要調整以保證 p∈［0，1］

for i in 1:n{

l［i］ ＜ -＊＊＊指定對數似然

p［i］＜ -l［i］-c 指定負二項分布參數 p 為 l［i］+c

y［i］ ＜-0對所有y設定為0

y［i］ ～ dnegbin(p［i］，1)}設定負二項分布似然

實例分析

麥克德里克(McKendrick)應用數學模型研究了1906年印度孟買一個村莊流行性霍亂傳染的數據〔4〕。數據顯示有大約75%的家庭沒有霍亂患者，有大約14%和7%的家庭有1和2個患者，患者人數超過3個的家庭數不足4%。

常數c的存在同樣不會影響貝葉斯后驗推斷。上式可看作Poisson分布數據yi=0，i=1，…，n的似然函數，參數λi=－li+c，c的作用在于保證λi＞0，可選一個大的正整數，如10000。WinBUGS程序如下:

c＜-10000設定常數c，在不合適時需要調整以保證λ＞0

for i in 1:n{

l［i］ ＜ -＊＊＊指定對數似然

lambda［i］＜ --l［i］+c 計算 Poisson 分布參數 λ

y［i］＜-0對所有y設定為0

y［i］～dpois(lambda［i］)}設定 Poisson 分布似然

(3)利用負二項分布指定非標準分布p表示家庭中成員不可能感染霍亂的概率，I()為指示函數，當x=0取值為1否則為0。顯然，ZIP的似然不能由WinBUGS軟件直接指定。ZIP的對數似然函數為

表2 1906年印度孟買村莊霍亂傳染數據的頻數分布

l［i］＜ -log(p*equals(x［i］，0)+(1-p)*exp(-mu)*pow(mu，x［i］)/exp(logfact(x［i］)))用ZIP的對數似然函數l［i］即可。WinBUGS函數equals(x，y)表示在 x=y時取值為1，否則為 0，pow(x，y)=xy，logfact(x)=log(x!)。Gibbs模擬總共抽樣15 000例，其中前5 000例作為 burn-in樣本拋棄不用〔6－7〕，則后驗樣本量為 10 000。p和 μ 的初始值分別設為0.5和1，c都指定為20。WinBUGS結果顯示三種方法的后驗統計量完全一樣，見表3。

表3 三種似然構建方法的ZIP模型后驗統計量

小 結

WinBUGS軟件基于Gibbs算法能夠從任意復雜的后驗分布模擬抽樣，極大地推廣了貝葉斯的應用，在WinBUGS軟件下能夠實現很多復雜的貝葉斯統計模型，如分層模型。本文在原有WinBUGS統計分布基礎上介紹了三種非標準分布的WinBUGS編程技巧，能夠滿足任意似然和先驗的分布類型，但都需要根據所采用分布的參數要求提供一個合適的常數c。c的作用在于保證所選用分布的參數符合實際意義，理論上c可以選擇一個相當大的數，在這種情況下一定能夠滿足分布的參數要求，但是在實際應用中，一個過大的c可能會導致計算機計算的浮點問題，這一點對二項和負二項分布尤其要注意，因此需要嘗試選擇一個滿足要求的數值同時保證不出現計算問題。例如ZIP模型中，二項和負二項分布中c選擇1000就會出現計算問題，但是對Poisson分布而言，只要c滿足要求其大小對計算并不會帶來什么實質問題，也基于這點原因，在實際應用時推薦使用Poisson分布構建任意分布的似然和先驗。

1．Spiegelhalter D，Thomas A，Best N，et al．WinBUGS user manual，Version 1．4，2003．http://www．mrc-bsu．cam．ac．uk/bugs/winbugs/manual14．pdf．

2．Gelman A，Carlin JB，Stern HS，et al．Bayesian data analysis，2nd ed．London:Chapman ＆ Hall，2004．

3．Ntzoufras I．Bayesian modeling using WinBUGS．New York:John Wiley ＆Sons，2009．

4．徐利治．現代數學手冊-隨機數學卷．武漢:華中科技大學出版社，1999．

5．Lambert D．Zero-inflated poisson regression with an application to defects in manufacturing．Technometrics，1992(34):1-14．

6．Gelfand AE，Smith AF．Sampling-based approaches to calculating marginal densities．Journal of the American Statistical Association，1990(85):398-409．

7．SASInstitute Inc．Preliminary capabilitites for bayesian analysis in SAS/STAT software．Cary，NC，USA，2006．