劉洋洋，陳 萍

(南京理工大學(xué) 理學(xué)院 南京 210094)

0 引 言

近年來，Bayes統(tǒng)計推斷取得了重大的發(fā)展，并且隨著Bayes方法的發(fā)展，該方法的優(yōu)點也得到了充分的體現(xiàn).通常Bayes方法的優(yōu)點可以概括為4點：它考慮了模型參數(shù)的先驗信息，并且參數(shù)的先驗信息用得越好，參數(shù)估計的精度就越高；該方法通常能夠較好地應(yīng)用于小樣本情形；能用抽樣的方法估計參數(shù)的后驗分布；能用抽樣的方法估計后驗分布的數(shù)字特征，如后驗期望、后驗方差、后驗眾數(shù)和后驗分位數(shù)等.在Bayes統(tǒng)計分析的問題中，蒙特卡洛抽樣(MCMC)方法的應(yīng)用比較廣泛，特別是MCMC方法中的M-H抽樣算法和Gibbs抽樣算法在參數(shù)估計的領(lǐng)域里運用得也比較多.本文利用這兩種算法針對廣義非線性模型的參數(shù)估計問題，提出從參數(shù)的條件后驗分布中抽取觀測值來估計參數(shù)值Bayes的方法，并用實例分析驗證了該方法的有效性和可行性.

1 廣義非線性模型

給定數(shù)據(jù)集(xi，yi)，i=1，2，…，n，設(shè)各yi之間相互獨立，E(yi)=μi，則廣義非線性模型可表示為

g(μi)=f(xi，β)

yi～ED(μi，σ2)，i=1，2，…，n

(1)

其中g(shù)(·)是嚴增可微函數(shù)，稱為聯(lián)系函數(shù)；ED(μi，σ2)表示指數(shù)族分布(陳希孺，1981；韋博成，2006)，yi=μ(xi，β)+εi，i=1，2，…，n，其中εi是均值為0的隨機誤差，即E(εi)=0.

yi的密度函數(shù)可表示為

p(yi；θi，φ)=

(2)

其中，θi稱為自然參數(shù)，σ2=φ-1稱為離差參數(shù)，并記μ=(μ1，μ2，…，μn)T，θ=(θ1，θ2，…，θn)T，根據(jù)指數(shù)族分布的性質(zhì)，有：

在廣義非線性模型中有廣泛應(yīng)用的偏差度概念，它相當于普通回歸模型中的殘差平方和，與尺度參數(shù)的估計有密切關(guān)系.對于廣義非線性模型式(1)，其偏差度函數(shù)定義為

d(yi，μi)=-2[yiθi-b(θi)-c(yi)]

因此式(1)可表示為β和σ2的函數(shù)，即

p(yi；μi(β)，σ2)=

(3)

由于各yi之間相互獨立，因此可得廣義非線性模型的概率密度函數(shù)為

2 模型的Bayes分析

2.1 分析方法

由Bayes公式可得模型參數(shù)的后驗密度為

p(β，σ2∣y，x)=

p(y∣β，σ2，x)p(β，σ2)=

(4)

由式(4)可知，要得到參數(shù)的后驗分布，要先指定參數(shù)的先驗分布.因為參數(shù)β是模型的回歸系數(shù)，σ2是隨機誤差分布的散度參數(shù)，所以通常假定p(β，σ2)=p(β)p(σ2)，假設(shè)p(β)和p(σ2)分別取為多元正態(tài)分布和逆Gamma分布，即假定參數(shù)β和σ2=φ-1的先驗分布為β～N[υ0，∑0]，φ～Gamma(α0，ρ0),故β和φ的聯(lián)合后驗分布可表示為

φα0-1exp{-φρ0}

因此，可得參數(shù)β的條件分布為

(β-υ0)T∑0-1(β-υ0)]}

(5)

參數(shù)φ的條件分布為

p(φ∣β，y，x)∝

(6)

如果a(yi，φ-1)為與φ無關(guān)的常數(shù)，則

如果a(yi，φ-1)∝φτexp {φ-k}，則

p(φ∣β，y，x)～Gamma(α0+nτ，nk+

(7)

由式(5)和式(6)可以看出，條件分布p(β∣φ，y，x)和p(φ∣β，y，x)一般都是一些非標準分布而且很復(fù)雜，不能直接從這些分布中抽取Bayes推斷所需要的觀測值，因此通常會采取MCMC方法中的Gibbs抽樣算法或M-H算法或這兩種算法相結(jié)合的混合算法來解決廣義非線性模型中參數(shù)估計的問題.

根據(jù)式(3)，Y=(y1，y2，…，yn)T的對數(shù)似然函數(shù)可表示為

由此可得

采用M-H算法從式(5)中抽取觀測樣本，設(shè)參數(shù)β的第j步迭代值為β(j)，則第(j+1)次迭代按如下步驟完成：

其中

采用M-H算法從式(6)中抽取觀測樣本，設(shè)φ的第j步迭代值為φ(j)，則第(j+1)次迭代按如下步驟完成：

若經(jīng)過計算發(fā)現(xiàn)有的參數(shù)條件分布可以用常見的標準分布表示出來，那么將選用Gibbs抽樣算法對其進行參數(shù)估計.Gibbs抽樣算法的一個優(yōu)勢就是所有抽樣數(shù)據(jù)由一維分布產(chǎn)生，樣本可以由標準統(tǒng)計分布生成；它不需要建議分布，每次接受新的樣本值.由此可知，Bayes估計的Gibbs抽樣算法步驟：

① 給定參數(shù)(β，φ)的初始值(β(0)，φ(0))，并令j=0；

② 已知φ(j)從后驗分布p(β∣φ(j)，y，x)中抽取觀測值β(j+1)；

③ 已知β(j+1)，從后驗分布p(φ∣β(j+1)，y，x)中抽取觀測值φ(j+1)；

④ 重復(fù)步驟②和步驟③得到β和φ的隨機樣本序列{β(j+1)，φ(j+1)：j=1，2，…，k}.

2.2 算法收斂性的診斷

無論使用哪一種抽樣方法，都要確定所得到的馬氏鏈的收斂性，即需要確定馬氏鏈達到收斂狀態(tài)時迭代的次數(shù)(達到收斂狀態(tài)前的那一段鏈稱為“預(yù)燒期”樣本).通常沒有一個全能的方法確定馬氏鏈的收斂性，監(jiān)視鏈的收斂性有許多方法，但是每種方法都是針對收斂性問題的不同方面提出的，因此，在絕大多數(shù)情況下，為了保證鏈的收斂性需應(yīng)用幾種不同的方法去診斷.常用的幾種診斷方法為蒙特卡洛誤差、樣本路徑圖、遍歷均值圖、自相關(guān)函數(shù)圖、Gelman-Rubin方法.本文采用的是利用樣本路徑圖、遍歷均值圖來判斷算法的收斂性.

樣本路徑圖是指用馬氏鏈迭代次數(shù)對生成的值作圖，若所有的值都在一個區(qū)域里且沒有明顯的周期性和趨勢性，那么可以假設(shè)收斂性已經(jīng)達到.

遍歷均值圖是一種很有用的圖方法，該方法是將馬氏鏈的累積均值對迭代次數(shù)作圖.這里累積均值是指此量直至當前迭代的平均值，如果累積均值在經(jīng)過一些迭代后基本穩(wěn)定，則表明算法已經(jīng)達到收斂.

2.3 Bayes推斷

本文采用的是M-H算法和Gibbs抽樣算法的混合算法，假設(shè)在第l次迭代時已經(jīng)收斂，則基于參數(shù)β和φ的觀測值序列{β(j)，φ(j)：j=1，2，…，k}，可得模型參數(shù)β和σ2=φ-1的Bayes估計.

(8)

(9)

令T=k-l，則由式(8)和式(9)給出的估計是相合估計，其相應(yīng)后驗協(xié)方差陣的相合估計為

因此，可以通過上述隨機樣本序列的樣本協(xié)方差陣的對角元素得到其對應(yīng)的標準誤差.

3 實例分析

通過實例分析，利用MCMC方法中Gibbs抽樣算法和M-H算法相結(jié)合的混合算法畫出參數(shù)的樣本路徑圖和均值遍歷圖，并求出廣義非線性模型中參數(shù)的Bayes估計，然后和已知的極大似然估計值進行比較，從而驗證該方法的有效性和可行性.

表1數(shù)據(jù)來自Whitmore(1986)，在此數(shù)據(jù)中，xi表示經(jīng)過市場調(diào)查得到的第i種產(chǎn)品的計劃銷售量，yi表示相應(yīng)的實際銷售量，Whitmore(1986)建議用以下逆高斯模型進行擬合：

(10)

表1 產(chǎn)品銷售數(shù)據(jù)

3.1 模型分析

在逆高斯模型式(10)中，結(jié)合式(2),有

φ=k=σ2

c(yi)=(2(yi))-1

s(yi，φ)=-{logφ-log(2πyi3)}

d(yi，μi)=(yi-(μi))2/μi2yi

假設(shè)參數(shù)β，γ，k服從的先驗分布分別為β～N(υ0，∑0)，γ～N(γ0，τ0)，k～Gamma(α0，ρ0).

根據(jù)式(5)可求得參數(shù)β和γ的條件分布，因為yi服從參數(shù)為μi和k=σ2的單形分布，故由式(7)可知，k的條件分布為Gamma分布，即

本例采用了混合抽樣的算法，即對參數(shù)β和γ采用M-H算法，對參數(shù)k采用Gibbs抽樣算法，然后用R軟件進行編程，對模型進行分析和求解.

3.2 參數(shù)的MCMC診斷

以下圖1和圖2為參數(shù)β，γ和k的樣本路徑圖和均值遍歷圖.

圖1 參數(shù)β和γ的MCMC診斷圖

圖2 參數(shù)k的MCMC診斷圖

由圖1中(a)(b)和圖2中的(e)可知參數(shù)β，γ和k的時間狀態(tài)圖顯示鏈混合得較好，波動比較穩(wěn)定，沒有明顯的周期性和趨勢性，這說明鏈達到了平穩(wěn)分布，由圖1中(c)(d)和圖2(f)可知參數(shù)β，γ和k的遍歷均值圖，從遍歷均值圖中可以看到去掉前8 000次迭代后所有的參數(shù)基本都趨于平穩(wěn)，這也表明了鏈的收斂性較好.因此丟掉前8 000次迭代值，選用之后的32 000次迭代值的平均值作為參數(shù)的Bayes估計，從圖2中可以看出鏈收斂后，參數(shù)均值的大致取值范圍.

3.3 參數(shù)估計

表2中為3個參數(shù)的Bayes估計值、標準差以及和極大似然估計值相比時的偏差.從偏差的值可以判定參數(shù)的Bayes估計和極大似然估計幾乎一致.

表2 產(chǎn)品銷售數(shù)據(jù)的Bayes估計

4 結(jié) 論

目前,研究廣義非線性模型參數(shù)估計問題的方法有很多種，用得較多的是似然估計方法和經(jīng)驗似然估計方法，每種方法都有它們各自的優(yōu)點.此處則選用了貝葉斯統(tǒng)計分析中蒙特卡洛抽樣方法中的M-H抽樣算法和Gibbs抽樣算法相結(jié)合的混合算法進行分析，求出參數(shù)的條件后驗分布后，根據(jù)其分布形式選出合適的算法，抽取每次迭代時的參數(shù)值，并利用參數(shù)的樣本路徑圖和均值遍歷圖驗證迭代時馬爾科夫鏈的收斂性；收斂后計算參數(shù)的后驗均值作為估計值，最后通過實證分析比較Bayes估計和極大似然估計的偏差，從而驗證M-H抽樣算法和Gibbs抽樣算法的簡潔性、有效性以及可行性.