徐登可 ，田瑞琴，吳劉倉

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

0 引言

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

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

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

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

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

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

(1)

其中：Yi表示第i個個體的實值響應變量，Zi表示p維標量型解釋變量向量，Xi(t)∈L2(T)是函數(shù)型解釋變量，L2(T)表示定義在概率空間上均值為零，二階矩有限的隨機過程，并且不失一般性假設T=[0,1]；θ=(θ1,θ2,…,θp)T是未知p維回歸參數(shù)，β(t)是定義在[0,1]上的平方可積函數(shù)；εi是獨立同分布服從于均值為零，方差為σ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)

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

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

(3)

(4)

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

(5)

1.2 似然函數(shù)

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

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

(6)

(7)

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

(8)

(9)

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

(10)

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

2 貝葉斯估計

2.1 先驗分布

為了應用貝葉斯估計方法，首先需要給出未知參數(shù)的先驗分布，具體為θ～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)合先驗分布為：

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

(11)

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

2.2 貝葉斯后驗推斷

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

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

(12)

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

?θ的滿條件分布：

(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 模擬研究

這部分通過2個隨機模擬例子來說明所提出的雙重部分函數(shù)型回歸模型和貝葉斯估計方法的有效性.

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

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

(17)

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

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

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

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

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

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

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

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

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

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

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

(18)

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

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

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

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

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

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

圖3 人均GDP的散點圖和對數(shù)的散點圖

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

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

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

5 結 論

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