非標準分布貝葉斯分析的WinBUGS實現(xiàn)*

徐州醫(yī)學院流行病與衛(wèi)生統(tǒng)計教研室(221002) 曾 平 王 婷 何 鵬

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

原理和方法

1．一個正態(tài)數(shù)據(jù)的WinBUGS程序

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

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

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

for i in 1:n{

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

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

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

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

}

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

表1 WinBUGS中的常用分布

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

c＜-10000設(shè)定常數(shù)c，在不合適時需要調(diào)整以保證 p∈［0，1］

for i in 1:n{

l［i］＜ -＊＊＊指定對數(shù)似然

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

y［i］＜-1對所有y設(shè)定為1

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

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

c＜-10000設(shè)定常數(shù)c，在不合適時需要調(diào)整以保證 p∈［0，1］

for i in 1:n{

l［i］ ＜ -＊＊＊指定對數(shù)似然

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

y［i］ ＜-0對所有y設(shè)定為0

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

實例分析

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

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

c＜-10000設(shè)定常數(shù)c，在不合適時需要調(diào)整以保證λ＞0

for i in 1:n{

l［i］ ＜ -＊＊＊指定對數(shù)似然

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

y［i］＜-0對所有y設(shè)定為0

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

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

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

l［i］＜ -log(p*equals(x［i］，0)+(1-p)*exp(-mu)*pow(mu，x［i］)/exp(logfact(x［i］)))用ZIP的對數(shù)似然函數(shù)l［i］即可。WinBUGS函數(shù)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和 μ 的初始值分別設(shè)為0.5和1，c都指定為20。WinBUGS結(jié)果顯示三種方法的后驗統(tǒng)計量完全一樣，見表3。

表3 三種似然構(gòu)建方法的ZIP模型后驗統(tǒng)計量

小 結(jié)

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

1．Spiegelhalter D，Thomas A，Best N，et al．WinBUGS user manual，Version 1．4，2003．http://www．mrc-bsu．cam．a(chǎn)c．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．徐利治．現(xiàn)代數(shù)學手冊-隨機數(shù)學卷．武漢:華中科技大學出版社，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．