M-H算法的建議分布選擇與比較

王丙參,魏艷華

(天水師范學院數學與統計學院,甘肅天水741001)

M-H算法的建議分布選擇與比較

王丙參,魏艷華

(天水師范學院數學與統計學院,甘肅天水741001)

文章比較研究了M一H算法與舍選法的聯系,給出了建議分布的選擇標準與幾種建議選擇,并比較了這幾種建議分布對馬氏鏈的影響,舉例說明了M一H算法在貝葉斯推斷中可行、穩定、有效。

M一H算法;舍選法;隨機數;建議分布;貝葉斯估計

0 引言

隨機模擬是研究復雜系統的常用手段,其難點在于目標分布的計算機仿真,即產生目標分布隨機數[1]。比如,產生多變量非標準形式且各分量之間相依分布的隨機數往往非常困難,難于實現,另外,在不完全數據場合進行統計推斷,往往需要產生給定截斷分布隨機數,但截斷分布參數的滿條件分布往往沒有顯示表達式,直接抽樣也很難。MCMC方法就是一種解決上述問題的簡單有效的貝葉斯計算方法[2,3],而M-H算法是MCMC的核心,已被列為影響工程技術與科學發展的十大算法之首,對其研究具有重要的實際意義[4]。在M-H算法具體操作中,建議分布選擇是否恰當直接關系到馬氏鏈的性質,比如混合性、收斂性等,進而影響后面各種統計計算[5,6]?；诖?本文比較研究了M-H算法與舍選法的聯系,給出了建議分布的選擇標準與幾種建議選擇,并比較了這幾種建議分布對馬氏鏈的影響,舉例說明了M-H算法在貝葉斯推斷中可行、穩定、有效。

1 對照舍選抽樣法研究M-H算法的基本思想與步驟

當某目標密度函數π(x)可被計算,但不易直接抽樣時,利用MCMC方法可獲得π(x)的一個近似樣本,基本思想是[5]:通過構造一個非周期不可約的馬氏鏈{X(t)，t=1，2，…},使其平穩分布恰好為目標分布π(x),這樣便可生成近似服從π(x)的樣本,從而進行各種統計推斷。M-H算法是重要性抽樣的實現,也是一種非常通用的構造馬氏鏈方法,它的動機是推廣舍選抽樣法。在舍選抽樣中,有目標分布π(x)、建議概率密度函數q(x)和一個接受準則h(x，y)(概率函數),同樣本文也假設π(x)是全樣本空間Ω上的目標分布,需要在Ω上產生樣本馬氏鏈{x1，…，xt，…},使它的穩定分布恰好為目標分布π(x)。首先,本文令q(x，y)表示由狀態x轉移到狀態y的概率,也常記為q(y|x),稱為建議分布或提案分布,則任意選擇不可約轉移概率q(x，y)及函數0＜α(.，.)≤1,對任意組合(x，y)，x≠y,則:

形成一個轉移核。注意,建議分布也常常用符號g(x，y)表示。

在t=0時刻,取X(0)=x(0),其中x(0)從初始分布h(x)中抽樣且滿足π(x(0))＞0。但是,初始點的選擇會影響馬氏鏈的性質,因此為了更快達到預期效果,在具體操作時,x(0)一般可取初始估計值,比如矩估計等。接著,給定X(t)=x,首先由q(.|x)產生潛在轉移y,然后以概率α(x，y)接受y(X(t+1)=y),以概率1-α(x，y)仍停留在x (X(t+1)=x)。

假設當前狀態X(t)=x(t),M-H算法具體步驟如下:

(1)從建議分布q(.|x(t))生成候選值x*。

(4)令t=t+1,返回第(1)步。

注意,M-H比率R(x(t)，x*)總是有定義,這是因為π(x)＞0且q(x，y)＞0。每完成上述一個循環,則生成一個目標分布隨機數。在具體實現M-H算法時,可選取多個初始點來檢驗得到的輸出結果是否一致,這也可看作是與最優化算法的結合。顯然,在M-H算法中,建議分布的主要目的是產生狀態轉移,即為每個狀態構造一個鄰域,并選中鄰域中某個鄰居狀態。要使得M-H算法產生的隨機數序列是馬氏鏈,這就需要要求建議分布q(x，y)在整個樣本空間Ω上有定義,通常,這是很容易實現的。舍選抽樣對所有的x有q(.|x)=q(.),是最簡單的情況,但MCMC方法推廣了舍選抽樣法,使建議分布變成了條件概率密度函數,而M-H算法的神奇之處也正在于此。當然,為了保證馬氏鏈狀態產生遍歷性的轉移,必須要求在x的某個鄰域Nx里成立:

M-H算法接下來要做的是,將選中的狀態與當前狀態相比,按一定的概率接受其中之一作為隨機序列的下一狀態。另一個與舍選抽樣法不同的是,接受概率不但取決于下一步狀態x(t+1),而且和當前狀態x(t)有關。

利用因M-H算法產生的隨機數序列,保證了X(t+1)只依賴X(t),而和以前產生的隨機數無關,所以M-H算法構造的鏈滿足馬氏性,即為馬氏鏈,但它是否是非周期不可約的則取決于建議分布q(x，y)的選擇。所謂的不可約是指馬氏鏈從任意狀態出發都能以一定的正概率轉移到其他狀態,即它的轉移是遍歷的。不具有遍歷性意味著馬氏鏈的整個狀態空間被分割成幾個互不相通的子空間,這樣,它的子空間的狀態就不能迭代到其他子空間中去,進而也不可能存在平穩分布。如果是非周期不可約的,則M-H算法構造的馬氏鏈具有唯一的極限平穩分布,且可以證明M-H算法構造的馬氏鏈以π(x)為其平穩分布。事實上,M-H算法產生的馬氏鏈是可逆的,即:

這表示馬氏鏈穩定時,由y→x的量等于由x→y的量,這就導致π(x)為其平穩分布。

利用M-H算法產生樣本的主要目的是估計X～π(x)的某一函數的期望。隨著t的增大,由M-H馬氏鏈產生隨機變量的分布越來越近似該鏈的平穩分布,故利用這種方法可以估計很多感興趣的量,如期望Eπ(h(X))、方差varπ(h(X))及尾概率P(X＞a),其中a為很大的數字,由馬氏鏈的極限性質可知,這些量的估計都是強相合的。由于生成的隨機數可以重復,并且這些抽樣點出現的頻率可以修正目標密度與提案分布密度的差異,所以本文應該保留馬氏鏈中的重復值。

但在實際應用中,我們不知道什么時候馬氏鏈才收斂到平穩分布。在M-H算法中,只有在極限情況下才有X(t)～π(x),但是,對于任何迭代都是有限的,不可能得到精確的邊際分布,而初始點的選擇對馬氏鏈是有影響的。馬氏鏈還未進入穩定階段之前的狀態稱為暫態階段,也可形象稱為熱機階段。于是,為了提高估計精度,我們通常舍去馬氏鏈的前D個值,即刪除一些初始生成值,也就是所謂的預燒期,但是關于預燒期的選擇沒有固定的標準,要根據具體情況具體分析,可能很大,比如D=20000,也可能很小,比如D=200。

2 建議分布的選擇與比較

通常,我們要求建議分布q(x，y)對于給定的x容易產生q(x，y)的隨機數,這樣才能使得抽樣變得簡單方便,否則,這只是將生成目標分布π(x)的困難轉移到了生成建議分布的困難,并不能解決問題。另外,一個具有某種特定性質的建議分布可以增強M-H算法的效果,并可直接通過接受概率的大小來反映。關于接受概率,應注意:接受概率要合適,具體問題要具體分布,而不是越大越好,因為太大會導致較慢的收斂性。Gelman等建議,當參數為1維時,接受概率應略小于0.5,當參數大于5維時,接受概率最好為0.25左右[5-7],當然,這只是個人建議,不同的學者可能有不同的看法。

不同的形式建議分布q(x，y)產生了相應的M-H算法變形,下面給出幾點具體建議供參考。

(1)M選擇:如果對任意的x，y,建議分布都有q(x，y)=q(y，x),即建議分布滿足對稱性,這時,M-H算法就是Metropolis算法(M算法),此時,接受概率M算法是MCMC方法的最早起源,由N.Metropolis等人于1953年在研究原子與分子的隨機運動問題時引入的隨機模擬方法。

特別有,如果q(x，y)=q(|x-y|),則成為隨機移動M算法。對稱建議分布很常用,例如,當給定x,q(x，y)可取正態分布,均值為x,協方差陣為常數陣,一種比較簡單的選擇就是取正態分布N(x(t)，σ2),關于σ的選擇,應注意:由大樣本性質,通常后驗分布都具有比較好的正態性,因此選擇正態分布作為建議分布是合適的,其均值應為上個狀態的值,而方差的大小決定了馬氏鏈在參數支撐集上的混合程度,因此σ的選擇決定了建議分布的好壞。當增量方差σ太大,大部分候選值被拒絕,這就導致M-H算法的效率太低,而σ太小,大部分候選值被接受,這時得到的馬氏鏈幾乎就是隨機游動,由于馬氏鏈從一個狀態到另一個狀態跨度太小,從而無法實現在整個支撐集上的快速移動,進而導致M-H算法的效率不高?？梢姦姨蠡蛱《紝е翸-H算法的效率較低。通常的做法是,在具體抽樣時監視候選值的接受概率,Robert,Gelman等建議最好控制在[0.15，0.5]。

(2)獨立抽樣:如果建議分布q(x，y)與當前狀態x無關,即q(x，y)=q(y),則由此建議分布導出的M-H算法稱為獨立抽樣。此時,由建議分布產生一個獨立鏈,抽取的每個候選值都與當前值獨立,M-H比率為并且如果g(x)＞0,π(x)＞0,則得到的馬氏鏈就是非周期不可約的。此處,M-H比率也可以表示重要比率,令其中π(x)是目標分布,q(x)為包絡分布,則這表明w(x(t))遠大于w(x*)的值時,馬氏鏈將會長時間停留在當前值上。進一步,我們有接受概率因此建議分布q(x)應與目標分布π(x)相似,且尾部包含π(x),這樣獨立抽樣的效果才可能好,這也導致獨立抽樣在實際中很少單獨使用。

(3)單變量M-H算法:如果目標分布X是高維的,則直接產生整個X比較困難,這時可以逐個抽取各分量,但是這需要用到完全條件分布。令…，Xn),考慮Xi|X～i,i=1，2，…，n的條件分布,選擇適當的轉移核q(xi→yi|x～i),對于給定的X～i=x～i產生下一個可能的yi,然后以概率:

(4)隨機游動鏈是通過簡單變化M-H算法得到的另一種馬氏鏈,實現方法如下:首先生成ε～h(ε),其中h為密度函數,令x*=x(t)+ε,就可得一個隨機游動鏈。在這種情況下,初始分布g(x*|x(t))=h(x*-x(t))。對于h,一般應選擇為正態分布N(0，σ2)、以原點為球心的球面上的均勻分布、尺度變換后的t分布。正態分布的方差可以控制候選值的擴散程度,如果目標分布是雙峰分布時,建議選擇σ可適當大點,在一般情況下,可令σ=1,即h～N(0，1)可以證明,如果h(x)在0的領域內為正且π(x)的支撐集連通,則生成鏈是非周期不可約的。

如果建議分布不隨時間t而改變,則稱M-H算法是齊次的。當然,本文也可構造非齊次的M-H算法,擊跑算法就是非齊次的M-H算法的典型代表,在這種方法中,從x(t)出發選擇的建議分布由兩步產生:先選擇移動方向,然后在此方向上產生移動距離。當X的狀態空間非常受限,且其他方法很難尋找所有空間的區域時,利用擊跑算法可能會取得比較好的效果。遺憾的是,非齊次的M-H算法收斂性質更難確定,甚至無法確定。

3 混合參數貝葉斯估計中M-H算法建議分布的比較

假定y為已知的觀測數據,待估參數θ的先驗分布為p(θ),則θ的貝葉斯估計往往依賴后驗分布p(θ|y)=cL (θ|y)p(θ),其中常數c為歸一化常數,未知,L(θ|y)為似然函數。由于在后驗分布p(θ|y)中c未知,所以它不能直接用于統計推斷,但是本文可利用中M-H算法獲得p(θ|y)的近似樣本,進而進行各種統計推斷。在M-H算法中,常常利用先驗分布作為建議分布,即建議分布g(θ)=p(θ),并取MH比率由于先驗分布的支撐集覆蓋后驗分布的支撐集,因此獨立鏈的平穩分布π(θ)即為后驗分布p(θ|y)。當然,先驗分布的選擇會顯著影響馬氏鏈的表現,如果選擇不合適,則會導致馬氏鏈收斂很慢,達不到預期效果。

假定觀測數據y1，…，y100來自混合正態分布:

觀測數據的直方圖見圖1,θ的先驗分布為U(0，1),本文可利用M-H算法構造平穩分布等于θ的后驗分布的馬氏鏈。已知觀測數據是由θ=0.7利用合成法生成的,故θ后驗密度應該集中在0.7附近,下面考察建議分布對馬氏鏈的影響。

圖1 由混合分布θN(6，0.52)+(1-θ)N(9，0.52)生成的100個觀測值直方圖

(1)利用Beta(1，1)(即U(0，1))作為建議分布,模擬結果如圖2。

圖2 不同初始值θ(1)生成的馬氏鏈軌跡(左圖θ(1)=0,右圖θ(1)=1)

圖2 表明馬氏鏈很快離開初始值,且似乎很容易從θ后驗支撐集的各個部分抽取值,這種表現稱為混合性良好。為了減少初始值的影響,本文省略前100次迭代,得到:左圖均值為0.6467,右圖均值為0.6483。可見,無論初始值為0,還是1,馬氏鏈都收斂,且樣本均值都與真值0.7非常接近,即模擬效果都不錯。

(2)利用Beta(2，10)作為建議分布,模擬結果如圖3。

圖3 不同初始值θ(1)生成的馬氏鏈軌跡(左圖θ(1)=0,右圖θ(1)=1)

圖3 生成的馬氏鏈慢慢離開初始值,由于擺動明顯,顯然,此鏈沒有收斂到其平穩分布,只能得到難以令人信服的結果,但是,其后驗分布仍是此鏈的極限分布,故長期運行此鏈,在原則上是可以估計θ的后驗分布。本文同樣省略前100次迭代,得到:左圖均值為0.5910,右圖均值為0.5966。可見,無論初始值為0,還是1,樣本均值都與真值0.7的差異比較大。

圖2和圖3說明建議分布對馬氏鏈的性質影響顯著,因此M-H算法的重點在于選擇合適的建議分布。

下面考慮用隨機游動鏈生成后驗分布。假定候選值θ*=θ(t)+ε,ε～U(-a，a),顯然,在馬氏鏈的迭代過程中,候選值有可能落在區間[0，1]外,這就導致迭代過程變得沒有實際意義。在這種情況下,本以為迭代結果很差,但是總體而言,迭代結果還是出乎意料的好。在實際操作中,馬氏鏈是如此的神奇,哪怕候選值不在區間[0，1],甚至密度函數值變為負值,迭代結果還是如此完美,仍在0.7附近波動,比如,當a=1時,迭代結果如圖4左圖所示,本文嘗試了很多取值,結果都是如此完美。但是,當a=8時,迭代結果就比較差了,如圖4右圖所示,這時候選值會經常不在[0，1]區間,但馬氏鏈基本不會接受不合理的候選值,因為不合理的候選值會導致的極大似然函數值變小。

圖4 θ*=θ(t)+ε且ε～U(-a，a)生成的馬氏鏈,左圖a=1,右圖a=8,初始值為0

解決候選值有可能落在區間[0，1]外的一種比較粗糙的解決辦法是,當θ?[0，1]時,重新令后驗值為0.5。一個更好的辦法是重新參數化。令現在,關于U生成一個隨機游動鏈,候選值u*=u(t)+ε1, ε1～U(-b，b)。下面給出兩種看待重新參數化的觀點。

(1)在θ空間運行馬氏鏈。這時,建議分布g(.|u(t))要通過變換成為θ空間的建議分布。令|J(θ(t))|表示θ→u變換的Jacobi行列式的絕對值在θ(t)的估值,則候選值θ*的M-H比率為:

(2)在u空間運行馬氏鏈。這時,θ的目標密度要變成u的目標密度,即則候選值U*=u*的M-H比率為:

顯然,在圖5(左)中,馬氏鏈混合性非常好,但是圖5 (右)中,效果就比較差,收斂速度較慢。為了減少初始值的影響,本文省略前100次迭代,得到:左圖均值為0.6476,標準差為0.0492,右圖均值為0.6502,標準差為0.0435,樣本均值都與真值0.7非常接近,右圖波動比左圖小,注意,這并不代表右圖的模擬結果更好,只是表明馬氏鏈混合性不好。

圖5 u空間運行馬氏鏈,左圖b=1,右圖b=10,初始值為0

圖6 u空間運行馬氏鏈,b=0.01,左圖初始值為0,右圖初始值為1

顯然,當b=0.01時,不管初始值為0,還是初始值為1,馬氏鏈運行效果都比較差。左圖的均值為0.5987,標準差0.0694,均值明顯偏離真實值。右圖的均值為0.7083,標準差0.0371,讀者也許很驚訝,均值竟然更接近0.7,認為模擬結果更好,其實不然,我們不能僅僅從一個偶然的好結果就判定整個方法好,還要考慮結果的產生過程,顯然,右圖馬氏鏈的收斂性和混合性都很差,模擬效果很差。

在所有方法中,估計值都在0.65附近,比真值0.7偏小,這主要是由于100個特定樣本值的性質導致的。

在重新參數化空間生成的隨機游動鏈與原始空間的生成的隨機游動鏈相比,具有很多不一樣的性質,重新參數化可以提高馬氏鏈的表現。

[1]王丙參,魏艷華,孫永輝.利用舍選抽樣法生成隨機數[J].重慶師范大學學報:自然科學版,2013,(6).

[2]魏艷華,王丙參,孫永輝.分組數據場合逆威布爾分布參數貝葉斯估計的混合Gibbs算法[J].四川大學學報:(自然科學版),2014,(1).

[3]魏艷華,王丙參.基于混合Gibbs算法定數截尾逆威布爾分布參數貝葉斯估計[J].統計與決策,2014,(11).

[4]肖柳青,周石鵬.隨機模擬方法與應用[M].北京:北京大學出版社, 2014.

[5]Givens GH,Hoeting JA.計算統計[M].王兆軍,劉民千,鄒長亮等譯.北京:人民郵電出版社,2009.

[6]Robert C P,Casella G.Monte Carlo Statistical Methods[M].北京:世界圖書出版社,2009.

[7]劉軍.科學計算中的蒙特卡羅策略[M].唐年勝,周勇,徐亮譯.北京大學出版社,2009.

(責任編輯/亦民)

O212

A

1002-6487(2016)20-0022-04

國家自然科學基金資助項目(61104045);天水師范學院中青年教師科研資助項目(TSA1404)

王丙參(1983—),男,河南南陽人,碩士,講師,研究方向:金融數學、統計計算。