雙重部分函數(shù)型回歸模型的貝葉斯估計(jì)

徐登可 ，田瑞琴，吳劉倉(cāng)

(1.杭州電子科技大學(xué) 經(jīng)濟(jì)學(xué)院, 浙江 杭州 310018; 2.杭州師范大學(xué) 數(shù)學(xué)學(xué)院, 浙江 杭州 311121; 3.昆明理工大學(xué) 理學(xué)院, 云南 昆明 650093)

0 引言

在現(xiàn)如今大數(shù)據(jù)時(shí)代，隨著電子信息技術(shù)的快速發(fā)展和先進(jìn)測(cè)量工具的出現(xiàn)，大量數(shù)據(jù)可以廉價(jià)地被收集和存儲(chǔ).其中在金融工程、環(huán)境科學(xué)、醫(yī)學(xué)、腦成像、公共衛(wèi)生等應(yīng)用領(lǐng)域，常常會(huì)獲得帶有明顯函數(shù)特性的數(shù)據(jù)，即觀測(cè)數(shù)據(jù)是在空間或時(shí)間上的一個(gè)或多個(gè)維度上獲得的，這一類型的數(shù)據(jù)稱之為函數(shù)型數(shù)據(jù).近年來(lái)，函數(shù)型數(shù)據(jù)分析已經(jīng)成為越來(lái)越熱的統(tǒng)計(jì)研究方向,且受到了很多統(tǒng)計(jì)學(xué)家的關(guān)注.例如，Ramsay和Dalzell[1]詳細(xì)介紹了具有函數(shù)型協(xié)變量和標(biāo)量響應(yīng)變量的函數(shù)型線性回歸模型.Shin[2]提出了部分函數(shù)型線性回歸模型，并研究了模型中未知回歸系數(shù)的理論性質(zhì).Lu等[3]研究了分位數(shù)部分函數(shù)型線性回歸模型，且獲得了模型中未知參數(shù)和函數(shù)型系數(shù)的漸近理論性質(zhì).Yu等[4]研究了變系數(shù)部分函數(shù)型線性分位數(shù)回歸模型.Zhou和Peng[5]研究了缺失數(shù)據(jù)下部分函數(shù)型線性回歸模型的參數(shù)估計(jì). Yu等[6]提出了單指標(biāo)部分函數(shù)型線性回歸模型.其他的函數(shù)型回歸模型研究還可以參見(jiàn)文獻(xiàn)[7-9].不難發(fā)現(xiàn)，上述函數(shù)型數(shù)據(jù)分析文獻(xiàn)主要基于模型方差齊性的假設(shè)，即模型誤差的方差是相等的.眾所周知，模型或者數(shù)據(jù)存在異方差時(shí)，采用現(xiàn)有的大多數(shù)統(tǒng)計(jì)推斷方法都有可能導(dǎo)致錯(cuò)誤的推斷.

然而，目前處理異方差數(shù)據(jù)最常用的方法是方差建模法，即不僅對(duì)均值建立回歸模型，同時(shí)也對(duì)方差建立回歸模型進(jìn)行分析，有些文獻(xiàn)稱之為雙重回歸模型或者聯(lián)合均值與方差模型.這個(gè)模型主要體現(xiàn)了對(duì)方差的重視，它能更好地解釋數(shù)據(jù)變化的原因和規(guī)律.特別最近這些年已經(jīng)有很多學(xué)者對(duì)基于方差建模的異方差模型研究了模型的參數(shù)估計(jì)、變量選擇以及異方差檢驗(yàn)等統(tǒng)計(jì)推斷.例如，吳劉倉(cāng)等[10]對(duì)聯(lián)合均值與方差模型提出一種同時(shí)對(duì)均值模型和方差模型的變量選擇方法；趙遠(yuǎn)英等[11]對(duì)響應(yīng)變量帶有不可忽略缺失數(shù)據(jù)的聯(lián)合均值與方差模型的貝葉斯估計(jì)問(wèn)題進(jìn)行了研究；戴琳等[12]基于聯(lián)合均值與方差模型研究了模型的參數(shù)估計(jì)與基于數(shù)據(jù)刪除模型考慮了統(tǒng)計(jì)診斷問(wèn)題.其他類似的相關(guān)研究還可以具體參見(jiàn)文獻(xiàn)[13-15].發(fā)現(xiàn)這些文獻(xiàn)大多數(shù)都是基于非函數(shù)型異方差模型展開(kāi)統(tǒng)計(jì)分析的，很少有文獻(xiàn)和學(xué)者基于方差建模研究異方差函數(shù)型數(shù)據(jù)回歸模型的貝葉斯參數(shù)估計(jì)等統(tǒng)計(jì)推斷問(wèn)題.

因此，本文主要基于方差建模針對(duì)異方差函數(shù)型數(shù)據(jù)提出了雙重部分函數(shù)型回歸模型，應(yīng)用Gibbs抽樣和Metropolis-Hastings算法相結(jié)合的混合MCMC算法研究模型的貝葉斯估計(jì).

1 模型與似然函數(shù)

1.1 雙重部分函數(shù)型回歸模型

考慮經(jīng)典部分函數(shù)型回歸模型如下：

(1)

其中：Yi表示第i個(gè)個(gè)體的實(shí)值響應(yīng)變量，Zi表示p維標(biāo)量型解釋變量向量，Xi(t)∈L2(T)是函數(shù)型解釋變量，L2(T)表示定義在概率空間上均值為零，二階矩有限的隨機(jī)過(guò)程，并且不失一般性假設(shè)T=[0,1]；θ=(θ1,θ2,…,θp)T是未知p維回歸參數(shù)，β(t)是定義在[0,1]上的平方可積函數(shù)；εi是獨(dú)立同分布服從于均值為零，方差為σ2的正態(tài)分布，即εi～N(0,σ2).

為方便描述，使用矩陣形式表示模型(1).令Y=(Y1,Y2,…,Yn)T，Z=(Z1,Z2,…,Zn)T，ε=(ε1,ε2,…,εn)T，X(t)=(X1(t),X2(t),…,Xn(t))T.模型(1)可以表示為：

(2)

其中：隨機(jī)誤差ε～N(0,σ2In),In是n維單位矩陣.

類似于Xu和Zhang[16],假設(shè)模型方差為異方差,并且利用一些解釋變量將方差參數(shù)建模為如下形式，即：

(3)

(4)

為了敘述方便，模型(4)可以簡(jiǎn)化為如下形式，包括均值模型和方差模型：

(5)

1.2 似然函數(shù)

令{(Yi,Zi,Hi,Xi),i=1,…,n}是來(lái)自模型(4)的獨(dú)立同分布樣本.定義函數(shù)型變量X(t)的協(xié)方差函數(shù)和經(jīng)驗(yàn)協(xié)方差函數(shù)分別為：

根據(jù)Karhunen-Loève表示定理可得：

(6)

(7)

因此，雙重部分函數(shù)型回歸模型可以近似為：

(8)

(9)

那么根據(jù)模型(9)可以得到似然函數(shù)：

(10)

其中V=Y-Zθ-Uγ.

2 貝葉斯估計(jì)

2.1 先驗(yàn)分布

為了應(yīng)用貝葉斯估計(jì)方法，首先需要給出未知參數(shù)的先驗(yàn)分布，具體為θ～N(θ0,Bθ)，γ～N(γ0,Bγ),δ～N(δ0,Bδ)，η～N(η0,Bη)，其中θ0,γ0,δ0,η0,Bθ,Bγ,Bδ和Bη是已知的超參數(shù).那么參數(shù)Θ=(θT,γT,δT,ηT)的聯(lián)合先驗(yàn)分布為：

π(θ,γ,δ,η)=p(θ)p(γ)p(δ)p(η)

(11)

其中p(θ)表示參數(shù)θ的先驗(yàn)概率密度函數(shù).

2.2 貝葉斯后驗(yàn)推斷

基于似然函數(shù)(10)和聯(lián)合先驗(yàn)分布(11)就可以獲得參數(shù)Θ=(θT,γT,δT,ηT)的聯(lián)合后驗(yàn)分布p(Θ|Y,X,H,Z)，具體如下：

p(Θ|Y,X,H,Z)∝L(θ,γ,δ,η|Y,Z,H,X)π(θ,γ,δ,η)

(12)

基于上式直接進(jìn)行抽樣和后驗(yàn)推斷是比較困難的.為了解決這個(gè)問(wèn)題，首先需要推導(dǎo)獲得每一個(gè)未知參數(shù)的滿條件分布，然后利用Gibbs抽樣和Metropolis-Hastings抽樣算法相結(jié)合的混合MCMC抽樣算法來(lái)從各自的滿條件分布中抽樣，具體如下.

?θ的滿條件分布：

(13)

?γ的滿條件分布：

(14)

?δ的滿條件分布：

p(δ|Y,Z,H,X,θ,γ,η)∝

(15)

?η的滿條件分布：

p(η|Y,Z,H,X,θ,γ,δ)∝

(16)

表1 未知參數(shù)Θ=(θ,γ,δ,η)的MCMC抽樣算法Tab.1 An MCMC-based sampling algorithm for unknown parameters Θ=(θ,γ,δ,η)

3 模擬研究

這部分通過(guò)2個(gè)隨機(jī)模擬例子來(lái)說(shuō)明所提出的雙重部分函數(shù)型回歸模型和貝葉斯估計(jì)方法的有效性.

3.1 例1：雙重部分函數(shù)型回歸模型

數(shù)據(jù)從如下雙重部分函數(shù)型回歸模型中產(chǎn)生：

(17)

Case I:選取好的先驗(yàn)信息θ0=(1,-0.5,0.5)T,Bθ=0.25×I3,δ0=(1,-0.5,0.5)T,Bδ=0.25×I3.

Case II:選取無(wú)先驗(yàn)信息θ0=(0,0,0)T,Bθ=10×I3,δ0=(0,0,0)T,Bδ=10×I3.

Case III:選取不精確的先驗(yàn)信息θ0=3×(1,-0.5,0.5)T,Bθ=I3,δ0=3×(1,-0.5,0.5)T,Bδ=I3.

其他超參數(shù)設(shè)置為γ0=0m,Bγ=10×Im,η0=0m,Bη=10×Im,這也表示選取比較弱的先驗(yàn)信息，其中0p表示全是0的p維向量.在模擬中分別令樣本量n=200,n=400和對(duì)于每一種情形下重復(fù)計(jì)算100次.

在上面設(shè)置的各種模擬環(huán)境下,應(yīng)用Gibbs抽樣和 Metropolis-Hastings算法相結(jié)合的混合MCMC算法來(lái)計(jì)算未知參數(shù)和函數(shù)型系數(shù)的貝葉斯估計(jì).對(duì)于每次重復(fù)產(chǎn)生的每一次數(shù)據(jù)集,MCMC算法的收斂性可以通過(guò)EPSR值來(lái)檢驗(yàn)[18],并且在每次運(yùn)行中觀測(cè)得到在 3 000 次迭代以后EPSR值都小于1.2.因此在每次重復(fù)計(jì)算中丟掉前 3 000 次迭代以后再收集J=2 000 個(gè)樣本來(lái)產(chǎn)生貝葉斯估計(jì).參數(shù)貝葉斯估計(jì)的模擬結(jié)果概括在表2～表3中. 另外，為了測(cè)量函數(shù)型系數(shù)估計(jì)的好壞，選擇用如下定義的RASE來(lái)衡量精確度：

表2 例1中當(dāng)樣本量n=200和不同的先驗(yàn)分布情況下未知參數(shù)的貝葉斯估計(jì)結(jié)果

表3 例1中當(dāng)樣本量n=400和不同的先驗(yàn)分布情況下未知參數(shù)的貝葉斯估計(jì)結(jié)果

表4 例1中在不同的樣本量下和不同的先驗(yàn)分布情況下函數(shù)型參數(shù)β(t)和α(t)的RASE的平均值和標(biāo)準(zhǔn)差

在表2～表3中,“Bias”表示基于100次重復(fù)計(jì)算未知參數(shù)的貝葉斯估計(jì)和真值之間的偏差,“SD”表示未知參數(shù)貝葉斯估計(jì)的標(biāo)準(zhǔn)差.從表2～表4中可以得到以下結(jié)論：1)在參數(shù)估計(jì)的偏差Bias和SD值方面,不管何種情形下貝葉斯估計(jì)都相當(dāng)精確，并且隨著樣本量的增大，模型中參數(shù)部分和函數(shù)型系數(shù)部分的貝葉斯估計(jì)結(jié)果變得越來(lái)越好.2)在不同的先驗(yàn)分布下，貝葉斯估計(jì)結(jié)果表現(xiàn)得都差不多，這也說(shuō)明提出的貝葉斯估計(jì)方法對(duì)先驗(yàn)分布的選取不是特別敏感.3)均值模型中參數(shù)估計(jì)的結(jié)果比方差模型中的參數(shù)估計(jì)效果要好一些.4)隨著樣本量的增大，RASE值的平均估計(jì)和標(biāo)準(zhǔn)差都變得越來(lái)越小，這也表明函數(shù)型系數(shù)估計(jì)得越來(lái)越好.從圖1和圖2中也展示了估計(jì)出來(lái)的函數(shù)型系數(shù)的曲線與相應(yīng)的真實(shí)函數(shù)的曲線逼近得都比較好，這與表4展示出來(lái)的結(jié)果是一樣的.總之,所有以上的模擬結(jié)果可以反映出所提出的貝葉斯估計(jì)方法能很好地恢復(fù)雙重部分函數(shù)型回歸模型中的真實(shí)信息.

3.2 例2：純粹的雙重函數(shù)型回歸模型

在這個(gè)例子中，從如下純粹的雙重函數(shù)型回歸模型中產(chǎn)生數(shù)據(jù)：

(18)

其中有關(guān)參數(shù)的設(shè)置和模擬環(huán)境與例1中一樣.另外，在這個(gè)例子中僅考慮無(wú)先驗(yàn)信息以及所有的貝葉斯分析結(jié)果展示在表5中，其中為了節(jié)省空間和避免累贅，函數(shù)系數(shù)曲線估計(jì)圖在此省略.從表中的結(jié)果可以發(fā)現(xiàn)，和預(yù)期的一樣,在這種模型下應(yīng)用提出的貝葉斯分析方法和例1中結(jié)果相似，且獲得的結(jié)果也是令人滿意的.這也說(shuō)明提出的貝葉斯估計(jì)方法也適應(yīng)于純粹的雙重函數(shù)型回歸模型.

表5 例2中在不同的樣本量下函數(shù)型參數(shù)β(t)和α(t)的RASE的平均值和標(biāo)準(zhǔn)差

圖1 例1中當(dāng)n=200和Case II先驗(yàn)信息下函數(shù)型系數(shù)β(t)和α(t)的真實(shí)函數(shù)曲線和平均估計(jì)曲線

圖2 例1中當(dāng)n=400和Case II先驗(yàn)信息下函數(shù)型系數(shù)β(t)和α(t)的真實(shí)函數(shù)曲線和平均估計(jì)曲線

4 實(shí)際數(shù)據(jù)分析

這部分將提出的雙重函數(shù)型回歸模型應(yīng)用到Growth數(shù)據(jù)集.該數(shù)據(jù)集描述了125個(gè)國(guó)家1961年到1985年的宏觀經(jīng)濟(jì)數(shù)據(jù).記Y表示1985年的人均GDP的對(duì)數(shù)，X(t)表示1961年到1985年年度儲(chǔ)蓄率，由真實(shí)投資與真實(shí)GDP比值得到.首先基于1985年的人均GDP及其對(duì)數(shù)的數(shù)據(jù)給出散點(diǎn)圖，如圖3所示，通過(guò)散點(diǎn)圖可以發(fā)現(xiàn)數(shù)據(jù)存在較顯著的異方差性.另外，在這主要想研究歷史的儲(chǔ)蓄率對(duì)1985年的人均GDP的影響，因此考慮如下雙重函數(shù)型回歸模型：

圖3 人均GDP的散點(diǎn)圖和對(duì)數(shù)的散點(diǎn)圖

其中采用無(wú)信息先驗(yàn)以及截?cái)鄥?shù)選取和模擬研究中一樣.另外,為了測(cè)試算法的收斂性,畫(huà)出了所有未知參數(shù)的EPSR值的圖, 且列在圖4中.從圖4中也能看出 3 000 次迭代以后所有參數(shù)的EPSR值都小于1.2，且接近1, 這表示 3 000 次迭代以后算法收斂了.在這里收集 3 000 次以后的后驗(yàn)樣本計(jì)算貝葉斯估計(jì).這樣就可以獲得函數(shù)型系數(shù)β(t)和α(t)的曲線估計(jì),如圖5所示.從圖5中可以看出，均值模型中函數(shù)型系數(shù)估計(jì)的曲線和方差模型中的函數(shù)型系數(shù)估計(jì)曲線很相似，總體上都隨著時(shí)間的推后而增加，在1982年達(dá)到最大，之后略有下降.這也表明歷史的儲(chǔ)蓄率對(duì)1985年的人均GDP總體上影響是正向的關(guān)系.

圖4 實(shí)際數(shù)據(jù)分析中所有參數(shù)的EPSR值Fig.4 EPSR values of all parameters in real data analysis

圖5 實(shí)際例子中函數(shù)型系數(shù)β(t)和α(t)的真實(shí)函數(shù)曲線和平均估計(jì)曲線

5 結(jié) 論

針對(duì)異方差函數(shù)型數(shù)據(jù)，本文基于方差建模的思想提出了雙重部分函數(shù)型回歸模型，其中使用函數(shù)型協(xié)變量對(duì)方差參數(shù)進(jìn)行建模.另外，運(yùn)用Karhunen-Loève表示定理來(lái)逼近函數(shù)型系數(shù)，以及基于給定的先驗(yàn)分布可以獲得未知參數(shù)的聯(lián)合后驗(yàn)分布和各個(gè)參數(shù)的條件分布，然后應(yīng)用Gibbs抽樣和Metropolis-Hastings算法相結(jié)合的混合MCMC算法來(lái)同時(shí)獲得均值模型和方差模型中未知參數(shù)和函數(shù)型系數(shù)的貝葉斯估計(jì).模擬研究顯示：1)隨著樣本量的增大，模型中參數(shù)部分和函數(shù)型系數(shù)部分的貝葉斯估計(jì)結(jié)果都是越來(lái)越好；2)貝葉斯估計(jì)方法對(duì)先驗(yàn)分布的選取不是特別敏感；3)均值模型中參數(shù)估計(jì)的結(jié)果比方差模型中的參數(shù)估計(jì)效果要好一些.另外，將提出的雙重函數(shù)型回歸模型應(yīng)用到Growth數(shù)據(jù)集，研究了歷史的儲(chǔ)蓄率對(duì)1985年的人均GDP的影響.兩個(gè)隨機(jī)模擬研究例子和實(shí)際數(shù)據(jù)分析都表明所提出的雙重部分函數(shù)型回歸模型和貝葉斯估計(jì)方法是可行有效的.