劉樂(lè)平，高 磊，楊 娜

（天津財(cái)經(jīng)大學(xué) 統(tǒng)計(jì)學(xué)系，天津300222）

一、引 言

1763年12月23日，理查德·普萊斯（Richard Price）在倫敦皇家學(xué)會(huì)會(huì)議上宣讀了托馬斯·貝葉斯的遺世之作，從此貝葉斯統(tǒng)計(jì)學(xué)誕生了。在整個(gè)19世紀(jì)，貝葉斯統(tǒng)計(jì)倍受爭(zhēng)議和冷落，主要原因是不知道如何恰當(dāng)?shù)靥幚硐闰?yàn)概率。20世紀(jì)上半葉，一個(gè)與貝葉斯統(tǒng)計(jì)完全不同的理論——頻率統(tǒng)計(jì)學(xué)（Frequentist statistics）迅速發(fā)展，在統(tǒng)計(jì)學(xué)領(lǐng)域幾乎占主導(dǎo)地位。盡管如此，貝葉斯統(tǒng)計(jì)思維仍在意大利的布魯諾·菲尼蒂和英國(guó)的哈羅德·杰弗里斯等少數(shù)統(tǒng)計(jì)學(xué)家的努力下，得以傳承［1］50－51。

現(xiàn)代貝葉斯運(yùn)動(dòng)的復(fù)興，始于20世紀(jì)下半葉，領(lǐng)軍人物是美國(guó)的吉米·薩維奇（Jimmy Savage）和英國(guó)的丹尼斯·林德利（Dennis Lindley），由于貝葉斯后驗(yàn)分布在高維計(jì)算上的困難，貝葉斯統(tǒng)計(jì)推斷非常難以實(shí)現(xiàn)，貝葉斯方法的應(yīng)用受到了很大的限制［2］。隨著計(jì)算機(jī)信息技術(shù)的發(fā)展和貝葉斯方法的改進(jìn)，特別是MCMC方法的發(fā)展、WinBUGS軟件以及R語(yǔ)言的廣泛應(yīng)用，原來(lái)復(fù)雜的數(shù)值計(jì)算問(wèn)題如今變得簡(jiǎn)單便捷，參數(shù)后驗(yàn)分布的模擬也趨于方便，現(xiàn)代貝葉斯分析的理論和應(yīng)用又重新得到了迅速的發(fā)展［3］。

MCMC（Markov Chain Monte Carlo）方法的發(fā)展，對(duì)現(xiàn)代貝葉斯分析的復(fù)興起著至關(guān)重要的作用。筆者在回顧MCMC方法發(fā)展歷史的基礎(chǔ)上，給出MCMC方法直觀(guān)性的實(shí)例，并介紹R軟件中MCMC的相關(guān)程序包，以此紀(jì)念貝葉斯定理250周年，因?yàn)镸CMC方法和現(xiàn)代貝葉斯兩者的發(fā)展和興起，不僅給出了解決問(wèn)題的方法，重要的是改變了考慮問(wèn)題的思路。

二、MCMC方法的起源

MCMC方法起源于軍事需要，來(lái)源于物理問(wèn)題，理論基礎(chǔ)是概率統(tǒng)計(jì)，研究手段是計(jì)算機(jī)技術(shù)。MCMC方法產(chǎn)生于第二次世界大戰(zhàn)期間，在研制原子彈的“曼哈頓計(jì)劃”過(guò)程中，1953年由美國(guó)物理學(xué)家Metropolis等人提出，經(jīng)加拿大統(tǒng)計(jì)學(xué)教授Hastings與其博士生于1970年推廣和完善。

（一）蒙特卡洛方法（Monte Carlo，MC）

蒙特卡洛方法誕生于第二次世界大戰(zhàn)期間的美國(guó)軍事研究重要基地——新墨西哥州的洛斯阿拉莫斯國(guó)家實(shí)驗(yàn)室。20世紀(jì)40年代中后期，在研制原子彈的“曼哈頓計(jì)劃”過(guò)程中，Stanislaw Ulam在病床上為了解決一個(gè)棘手的組合計(jì)算問(wèn)題（紙牌游戲“solitaire”中獲勝概率的計(jì)算），產(chǎn)生了蒙特卡洛方法的最初理念。約翰·馮·諾伊曼（John von Neumann）將這種思路直接應(yīng)用到核問(wèn)題的中子擴(kuò)散研究中［4］。Nicholas Metropolis建議將這種方法命名為“蒙特卡洛”①之所以用賭城"蒙特卡洛"來(lái)命名這種隨機(jī)模擬算法，Metropolis（1987）說(shuō)是由于烏拉姆（Ulam）的叔叔特別沉迷于摩納哥的大賭場(chǎng)——Monte Carlo Casino（http：／／en．wikipedia．org／wiki／Monte＿Carlo＿method）。。

1946年2月，經(jīng)過(guò)三年的研發(fā)，世界上第一臺(tái)計(jì)算機(jī)——ENIAC誕生。ENIAC的問(wèn)世與蒙特卡洛方法的發(fā)展密切相關(guān)。1947年，約翰·馮·諾依曼對(duì)蒙特卡洛方法進(jìn)行了拓展，以此來(lái)解決熱核和裂變問(wèn)題。同年，Ulam和約翰·馮·諾依曼發(fā)明了逆向和接受－拒絕（accept－reject）技術(shù)來(lái)模擬均勻分布。1949年，在由蘭德、國(guó)家統(tǒng)計(jì)局和橡樹(shù)嶺實(shí)驗(yàn)室（Rand，NBS and the Oak Ridge laboratory）支持的關(guān)于蒙特卡洛的研討會(huì)上，Metropolis和Ulam共同發(fā)表了第一篇關(guān)于蒙特卡洛方法的論文［5］。

（二）MCMC方法的開(kāi)篇之作

MCMC算法的起源與世界上第二臺(tái)計(jì)算機(jī)——MANIAC有關(guān)。1952年初，在 Metropolis指導(dǎo)下，MANIAC在洛斯阿拉莫斯國(guó)家實(shí)驗(yàn)室誕生。Nicolas Metropolis是創(chuàng)建洛斯阿拉莫斯國(guó)家實(shí)驗(yàn)室的15名科學(xué)家之一（1915年生于美國(guó)，1941年獲得芝加哥大學(xué)物理學(xué)博士學(xué)位，1943年4月來(lái)到洛斯阿拉莫斯國(guó)家實(shí)驗(yàn)室，直到1999年去世，在實(shí)驗(yàn)室工作了56年），Nicolas Metropolis既是物理學(xué)家又是數(shù)學(xué)家，既是美國(guó)科學(xué)院的院士又是美國(guó)數(shù)學(xué)會(huì)和物理學(xué)會(huì)的會(huì)士（Fellow）。20世紀(jì)50年代初，借助改進(jìn)的計(jì)算機(jī)設(shè)備，Micolas Metropolis他與Ulam一起合作，負(fù)責(zé)設(shè)計(jì)研發(fā)氫彈［5］。

1953年6月，在化學(xué)物理期刊（Journal of Chemical Physics）上Metropolis等人發(fā)表了MCMC方法的開(kāi)篇之作《通過(guò)快速計(jì)算機(jī)器計(jì)算狀態(tài)方程》（Equations of state calculations by fast computing machines），論文的主要關(guān)注點(diǎn)是在R2N內(nèi)計(jì)算以下積分公式（類(lèi)似于貝葉斯的后驗(yàn)分布）②此文章到2013年9月10日被引用次數(shù)為25 191次。：

其中θ表示R2里的N個(gè)粒子，粒子的能量E表示為：

其中V是一個(gè)勢(shì)函數(shù)，dij是θ中粒子i和j之間的歐氏距離。

通過(guò)溫度T、玻爾茲曼（Boltzmann）常數(shù)k和正規(guī)因子

可以將波爾茲曼分布exp｛－E（θ）／kT｝參數(shù)化。

由于θ是2N維向量，所以不可能用數(shù)值積分計(jì)算。對(duì)于高維計(jì)算問(wèn)題，由于exp｛－E（θ）／kT｝對(duì)于粒子系統(tǒng)隨機(jī)配置的大多數(shù)實(shí)現(xiàn)值結(jié)果都非常小，所以標(biāo)準(zhǔn)的蒙特卡洛方法不能正確地逼近積分值犐。

為提高蒙特卡洛方法的有效性，Metropolis等人（1953）提出了N個(gè)粒子隨機(jī)游走的修正方法，即對(duì)N個(gè)粒子中的每個(gè)粒子（X，Y）都做一個(gè)微小的隨機(jī)移動(dòng)，設(shè)ξ1，ξ2是（－1，1）上的隨機(jī)數(shù)，對(duì)X和Y分別做一個(gè)隨機(jī)的變換，令：

產(chǎn)生新粒子（X′，Y′），并且假定新粒子以概率

被接受，否則將重復(fù)原粒子。新粒子被接受后，可以計(jì)算新粒子結(jié)構(gòu)和原粒子結(jié)構(gòu)的能量差ΔE。要特別注意的是，Metropoli等人文章中一次只移動(dòng)一個(gè)粒子，而不是一起移動(dòng)所有的粒子，因此使這種算法看起來(lái)像Gibbs樣本法的雛形［6］。

（三）Metropolis算法的推廣

Nicholas Metropolis是物理學(xué)家，他在核物理研究中碰到的粒子分布高緯計(jì)算問(wèn)題，對(duì)統(tǒng)計(jì)學(xué)提出了挑戰(zhàn)。Hastings（1970）和他指導(dǎo)的唯一博生生Peskun（1973）圓滿(mǎn)地解決了這個(gè)難題，攻克了常規(guī)蒙特卡洛方法遇到的維度瓶頸問(wèn)題，將Metropolis算法一般化，并推廣為一種統(tǒng)計(jì)模擬工具，形成了Metropolis－Hastings（M－H）算法。W．Keith Hastings1930年出生于加拿大，1962年博士畢業(yè)于多倫多大學(xué)數(shù)學(xué)系（當(dāng)時(shí)統(tǒng)計(jì)學(xué)在數(shù)學(xué)系中），1964年在貝爾實(shí)驗(yàn)室工作兩年后，在母校數(shù)學(xué)系任統(tǒng)計(jì)教職，他的研究興趣主要關(guān)注概率分布隨機(jī)抽樣的蒙特卡洛方法①關(guān)于Hastings的資料較少，較詳細(xì)的網(wǎng)絡(luò)資料見(jiàn)http：／／probability．ca／hastings／。。當(dāng)時(shí)多倫多大學(xué)化學(xué)系John Valleau教授的研究團(tuán)隊(duì)面臨一個(gè)“粒子系統(tǒng)的平均勢(shì)能估計(jì)難題”，即在一個(gè)確定勢(shì)場(chǎng)（defined potential field）中，100個(gè)粒子組成了一個(gè)粒子系統(tǒng)，由于每個(gè)粒子有6維坐標(biāo)，所以粒子系統(tǒng)就有600個(gè)維度。那么，如何使用Metropolis方法估計(jì)高緯粒子系統(tǒng)的平均勢(shì)能。當(dāng)Hastings得知這個(gè)問(wèn)題可以利用馬爾可夫鏈、并通過(guò)抽取高緯分布的隨機(jī)樣本輕松解決時(shí)，他意識(shí)到這個(gè)方法對(duì)于統(tǒng)計(jì)學(xué)的重要意義，于是便全身心投入此方法及其變化的研究中，完成了MCMC方法史上里程碑的論文——Monte Carlo Sampling Methods Using Markov Chains and their Applications（《馬爾可夫鏈蒙特卡洛抽樣方法及其應(yīng)用》），發(fā)表在1970年的Biometrika上［7］。

M－H算法相對(duì)于Metropolis方法而言，看起來(lái)更像是專(zhuān)業(yè)的統(tǒng)計(jì)模擬工具，最重要的是，M－H算法不要求“建議分布函數(shù)”必須是對(duì)稱(chēng)的，從而應(yīng)用起來(lái)更加靈活方便。Hastings在文章里舉了三個(gè)例子：第一個(gè)目標(biāo)分布是泊松分布，采用的建議分布是隨機(jī)游走；第二個(gè)目標(biāo)分布是正態(tài)分布，建議分布是均勻隨機(jī)游走，但此均勻分布的中心不是馬氏鏈的當(dāng)前值θt，而是－θt；第三個(gè)目標(biāo)分布是多元分布，Hastings引進(jìn)了Gibbs抽樣的策略，每次只更新其中一個(gè)變量，這樣的轉(zhuǎn)移矩陣同樣滿(mǎn)足平穩(wěn)條件，因?yàn)槊看味茧x開(kāi)了目標(biāo)不變變量。三年后，Peskun發(fā)表了題為Optimum Monte－Carlo Aampling Using Markov Chains（《最優(yōu)馬爾可夫鏈蒙特卡洛抽樣方法》）的文章，比較了Metropolis和Barker的接受概率的形式，也證明了在離散情形下Metropolis是最優(yōu)的選擇，這里的最優(yōu)性可以通過(guò)經(jīng)驗(yàn)平均值的漸近方差來(lái)表示［8］。

三、MCMC方法的發(fā)展

從Peskun之后，在統(tǒng)計(jì)學(xué)研究領(lǐng)域里關(guān)于MCMC方法的研究沉寂了近十年之久。隨后，在模式識(shí)別、圖像分析和空間統(tǒng)計(jì)學(xué)等領(lǐng)域中，出現(xiàn)了關(guān)于 MCMC方法應(yīng)用的文章，其中 Geman and Geman發(fā)表了在MCMC方法史上具有重要突破性的文章［9］。該文基于隨機(jī)松弛（Stochastic relaxation）算法，采用Gibbs分布對(duì)圖像的貝葉斯恢復(fù)進(jìn)行了研究，提出了Gibbs采樣的概念并將其引入到統(tǒng)計(jì)應(yīng)用領(lǐng)域，Robert和Casella（2011）將此文稱(chēng)為“革命的種子”，吹響了MCMC方法革命的號(hào)角。

1987年，Tanner和 Wong在論文中采用基于多個(gè)條件分布進(jìn)行模擬的方法，這種思路等價(jià)于從聯(lián)合分布進(jìn)行模擬，基本具備了Gibbs采樣的雛形。這篇文章被選入了美國(guó)統(tǒng)計(jì)學(xué)會(huì)雜志（Journal of the American Statistical Association）的 討 論 論文［10］。需要指出的是，Tanner和 Wong的論文雖然已經(jīng)基本具備了Gelfand和Smith（1990）文章的雛形，但與后者相比，其影響有限。原因之一是Tanner和Wong的論文似乎只是應(yīng)用到缺失數(shù)據(jù)問(wèn)題，特別是論文的題目中就有“數(shù)據(jù)補(bǔ)全”這樣的字眼；另一原因是Tanner和Wong論文的理論基礎(chǔ)不是馬爾可夫鏈，而是函數(shù)分析，特別是要求馬爾可夫轉(zhuǎn)移核一致有界而且連續(xù)，也許正因如此影響了很多數(shù)學(xué)背景不強(qiáng)的潛在研究人員。

一系列以MCMC方法和貝葉斯為主題的學(xué)術(shù)會(huì)議，展示了Gibbs抽樣爆炸式的發(fā)展過(guò)程。1986年夏，Adrian Smith做了關(guān)于分層模型（Hierarchical models）的系列學(xué)術(shù)演講；1989年6月，在魁北克省舍布魯克市（Sherbrooke，Québec）舉行的貝葉斯學(xué)術(shù)會(huì)議上，Adrian Smith第一次詳細(xì)闡釋了Gibbs抽樣的本質(zhì)，這種方法的廣度與深度震撼了與會(huì)者。1990年，Gelfand和Smith發(fā)表論文Sampling－based approaches to calculating marginal densities（《基于抽樣的邊際分布計(jì)算方法》），將這種思想闡述得更為深刻和完整，成為主流統(tǒng)計(jì)學(xué)界大規(guī)模使用MCMC方法的真正起點(diǎn)［11］。

20世紀(jì)90年代，是MCMC方法發(fā)展的黃金時(shí)期，并在理論研究方面獲得了很多突破：1991年，Alan Gelfand，Prenatal Goel和 Adrian Smith在俄亥俄州立大學(xué)（Ohio State University）舉辦了 MCMC方法會(huì)議，相關(guān)討論成果都已成為MCMC領(lǐng)域非常有影響力的論文。1992年的5月，皇家統(tǒng)計(jì)學(xué)會(huì)（Royal Statistical Society）召開(kāi)了關(guān)于“Gibbs抽樣與其他MCMC方法”的會(huì)議，有四篇論文得到了眾多學(xué)者的重視，并發(fā)表在JRSSB 1993年的第一期上。

在理論研究獲得飛速發(fā)展的同時(shí)，MCMC的應(yīng)用研究也取得了可喜的成果，MCMC方法之所以發(fā)展迅速，一個(gè)重要的原因就是它的簡(jiǎn)便實(shí)用。以Gelfand和Smith提出的非常簡(jiǎn)單的隨機(jī)效應(yīng)模型為例：

其中θi～N（μ，σθ2），εij～N（0，σε2），獨(dú)立于θi。

在頻率學(xué)派看來(lái)，對(duì)方差分量進(jìn)行估計(jì)比較困難；對(duì)于傳統(tǒng)貝葉斯學(xué)派來(lái)說(shuō)，這也是一個(gè)噩夢(mèng)，因?yàn)榇藭r(shí)的積分難以求得；對(duì)以MCMC為特征的現(xiàn)代貝葉斯分析來(lái)說(shuō)，在給定參數(shù) μ，σθ2，σε2先驗(yàn)信息的情況下，問(wèn)題可以很方便地通過(guò)Gibbs抽樣解決?；谠撾S機(jī)效應(yīng)模型，學(xué)者們利用Gibbs抽樣迅速得到了貝葉斯線(xiàn)性混合模型的參數(shù)估計(jì)，該隨機(jī)效應(yīng)模型的其他應(yīng)用還包括空間統(tǒng)計(jì)、變點(diǎn)分析、回歸中的變量選擇和基因組學(xué)等。

到了20世紀(jì)90年代中期，MCMC理論漸趨成熟，理論突破開(kāi)始變緩，但仍然取得了一些成果。許多動(dòng)態(tài)系統(tǒng)可以用狀態(tài)空間模型（state－space model）來(lái)描述，而序貫蒙特卡洛模擬（Sequential Monte Carlo，SMC）是分析狀態(tài)空間模型的重要工具。Roberts和 Rosenthal（2007）建立了SMC與MCMC的聯(lián)系，提出了適應(yīng)性 MCMC［12］。

Green（1995）提出的逆跳的馬爾可夫鏈蒙特卡洛模 擬 （Reversible Jump Markov Chain Monte Carlo，RJMCMC），可被視為MCMC方法第二次革命的開(kāi)端，所構(gòu)建的馬氏鏈不僅可以在一個(gè)模型的參數(shù)空間內(nèi)進(jìn)行轉(zhuǎn)移，還可以在不同模型（維度可以不同）之間跳躍，從而為貝葉斯模型選擇提供了強(qiáng)大工具［13］。Richardson和 Green（1997）將 RJMCMC應(yīng)用到混合階估計(jì)中；Brooks，Giudici和Roberts提出了提高RJMCMC轉(zhuǎn)移效率的方法；Marin和Robert將RJMCMC應(yīng)用到變量選擇中。由于MCMC方法得到的并非獨(dú)立樣本，具有自相關(guān)性，因此不能直接用于參數(shù)估計(jì)，尤其是不能用于直接估計(jì)參數(shù)的標(biāo)準(zhǔn)誤差。為了克服自相關(guān)性，Hobert等人提出在原始馬氏鏈轉(zhuǎn)移鏈基礎(chǔ)上間隔抽取構(gòu)成新的樣本序列，從而克服自相關(guān)性，可以得到參數(shù)估計(jì)以及估計(jì)誤差［14］。

四、MCMC方法直觀(guān)性實(shí)例

MCMC方法實(shí)質(zhì)上就是利用馬爾可夫鏈進(jìn)行蒙特卡洛模擬。它可以分解成兩個(gè) MC：前一個(gè)MC是馬爾可夫鏈 （Markov Chain），后一個(gè) MC是蒙特卡洛方法 （Monte Carlo）。在此，按照蒙特卡洛方法、馬爾可夫鏈和Metropol－Hastings算法的次序，介紹MCMC方法的直觀(guān)性實(shí)例。

（一）蒙特卡洛方法（Monte Carlo）：圓周率π的估計(jì)

蒙特卡洛方法與數(shù)值計(jì)算中的其他確定性算法不同，是一種以概率統(tǒng)計(jì)理論為基礎(chǔ)的非確定性隨機(jī)模擬算法，也稱(chēng)統(tǒng)計(jì)模擬方法。它使用隨機(jī)數(shù)（或偽隨機(jī)數(shù)）解決很多不確定性的計(jì)算問(wèn)題，將所求解的問(wèn)題同一定的概率分布相聯(lián)系，并用電子計(jì)算機(jī)實(shí)現(xiàn)統(tǒng)計(jì)模擬或抽樣，從而獲得問(wèn)題的近似解。一般來(lái)說(shuō)，蒙特卡洛方法可以粗略地分成兩類(lèi)：

一類(lèi)是所求解的問(wèn)題本身具有內(nèi)在的隨機(jī)性，借助計(jì)算機(jī)的運(yùn)算能力可以直接模擬這種隨機(jī)的過(guò)程。例如在核物理研究中，分析中子在反應(yīng)堆中的傳輸過(guò)程。中子與原子核作用受到量子力學(xué)規(guī)律的制約，人們只能知道它們相互作用發(fā)生的概率，卻無(wú)法準(zhǔn)確獲得中子與原子核作用時(shí)的位置以及裂變產(chǎn)生的新中子的行進(jìn)速率和方向?？茖W(xué)家依據(jù)其概率進(jìn)行隨機(jī)抽樣得到裂變位置、速度和方向，進(jìn)而在模擬大量中子的行為后經(jīng)過(guò)統(tǒng)計(jì)就能獲得中子傳輸?shù)姆秶源俗鳛榉磻?yīng)堆設(shè)計(jì)的依據(jù)。

另一類(lèi)是所求解的問(wèn)題可以轉(zhuǎn)化為某種隨機(jī)分布的數(shù)字特征值，比如隨機(jī)事件出現(xiàn)的概率，或者隨機(jī)變量的期望值。通過(guò)隨機(jī)抽樣的方法，以隨機(jī)事件出現(xiàn)的頻率估計(jì)其概率，或者以樣本的數(shù)字特征估算隨機(jī)變量的數(shù)字特征，并將其作為問(wèn)題的解。這種方法多用于求解復(fù)雜的多維積分問(wèn)題（http：／／zh．wikipedia．org／）。

例如，假設(shè)要計(jì)算一個(gè)不規(guī)則圖形的面積，圖形的不規(guī)則程度和分析性計(jì)算（比如積分）的復(fù)雜程度成正比，與利用定積分方法計(jì)算不規(guī)則圖形面積的方法（任意細(xì)分后求和再取極限得到精確值）不同，蒙特卡洛方法的簡(jiǎn)單直觀(guān)思路為：假想有一袋小玻璃珠子，把玻璃珠子均勻地朝這個(gè)圖形上撒，再數(shù)此圖形中有多少顆玻璃珠子，而玻璃珠子的數(shù)目就近似代表圖形的面積。當(dāng)玻璃珠子越小、撒得越多的時(shí)候，結(jié)果就越精確。借助計(jì)算機(jī)程序可以生成大量均勻分布坐標(biāo)點(diǎn)（代替玻璃珠子），然后統(tǒng)計(jì)出圖形內(nèi)的點(diǎn)數(shù)，通過(guò)它們占總點(diǎn)數(shù)的比例和坐標(biāo)點(diǎn)生成范圍的面積，就可以求出圖形面積。

利用蒙特卡洛方法，可以非常簡(jiǎn)便地近似計(jì)算圓周率π：在一個(gè)邊長(zhǎng)為單位1的正方形內(nèi)，畫(huà)出一個(gè)半徑為1的1／4圓，見(jiàn)圖1。

圖1 利用蒙特卡洛方法計(jì)算圓周率圖

向正方形內(nèi)隨機(jī)投點(diǎn)，設(shè)總投點(diǎn)數(shù)為N；統(tǒng)計(jì)圓內(nèi)的點(diǎn)數(shù)，用n表示，則n與總點(diǎn)數(shù)N的比值就近似等于1／4圓的面積（π／4）和正方形面積（1）之比，即π／4，所以π≈4n／N。

通過(guò)計(jì)算機(jī)來(lái)模擬上述隨機(jī)投點(diǎn)過(guò)程：首先，每次隨機(jī)生成兩個(gè)服從均勻分布的0到1之間的隨機(jī)數(shù)，設(shè)為x和y，分別以x和y為橫縱坐標(biāo)，就得到圖中的一個(gè)隨機(jī)投點(diǎn)；其次，再判斷這個(gè)點(diǎn)是否在圓內(nèi)，而固定x只需判別y是否小于或等于即可；最后，生成N個(gè)隨機(jī)點(diǎn)，加總圓內(nèi)點(diǎn)的總數(shù)n，計(jì)算4n／N即可。當(dāng)隨機(jī)點(diǎn)N取值越大時(shí)，其結(jié)果越接近于圓周率（當(dāng)投點(diǎn)達(dá)到3萬(wàn)次時(shí)，用蒙特卡洛方法可以得到π≈3．143 6）。

（二）馬爾可夫鏈 （Markov Chain）：股市的表現(xiàn)

理解了蒙特卡洛方法的基本思路，再看馬爾可夫鏈及其平穩(wěn)分布的相關(guān)理論。馬爾可夫鏈（以下簡(jiǎn)稱(chēng)馬氏鏈）的定義如下：

簡(jiǎn)言之，即馬氏鏈的下一個(gè)狀態(tài)只依賴(lài)于當(dāng)前狀態(tài)，而與已經(jīng)過(guò)去的狀態(tài)沒(méi)有關(guān)系。它是介于獨(dú)立和相關(guān)之間的一種理想形式，例如在研究時(shí)間序列時(shí)，為了便于分析，希望昨天、今天和明天表示的三種隨機(jī)狀態(tài)（或變量）完全獨(dú)立，但實(shí)際上可能它們完全相關(guān)，變量之間相互獨(dú)立太理想，完全相關(guān)又太復(fù)雜，馬氏鏈對(duì)此進(jìn)行了簡(jiǎn)化，假設(shè)明天只與今天相關(guān)，而與昨天無(wú)關(guān)（獨(dú)立）①馬氏鏈直觀(guān)的解釋是荷花池中青蛙的跳躍過(guò)程；有趣的解釋是，芭蕾舞團(tuán)在招考少兒芭蕾舞女演員時(shí)，預(yù)計(jì)她未來(lái)的體型，只需看她的母親而不需要看她的姥姥。。換言之，即在馬氏鏈過(guò)程中，在給定當(dāng)前知識(shí)或信息的情況下，只有當(dāng)前的狀態(tài)用來(lái)預(yù)測(cè)將來(lái)，過(guò)去（即當(dāng)前以前的歷史狀態(tài)）對(duì)于預(yù)測(cè)將來(lái)（即當(dāng)前以后的未來(lái)狀態(tài)）是無(wú)關(guān)的。

下面以來(lái)自wikipedia的例子來(lái)說(shuō)明馬氏鏈及其平穩(wěn)性：假設(shè)股票市場(chǎng)在一周內(nèi)的表現(xiàn)有三種可能狀態(tài)：熊市—狀態(tài)1、牛市—狀態(tài)2、平衡市—狀態(tài)3。股市行情的轉(zhuǎn)化見(jiàn)圖2。

圖2 馬氏鏈三種狀態(tài)轉(zhuǎn)移示意圖

由圖2可見(jiàn)，如果當(dāng)前股市為牛市，則下周為熊市的概率為7．5%、平衡市的概率為2．5%、依然為牛市的概率為90%，其他行情轉(zhuǎn)化與之類(lèi)似，則該馬氏鏈的一步轉(zhuǎn)移概率矩陣為：

考慮三種初始情況：當(dāng)前行情為牛市、當(dāng)前行情為熊市、當(dāng)前行情為平衡市，它們對(duì)應(yīng)的初始概率分布為：π0＝（1，0，0）、π0＝（0，1，0）、π0＝（0，0，1）。分析在這三種情況下股市行情在未來(lái)30個(gè)星期之內(nèi)的轉(zhuǎn)移變化情況，見(jiàn)圖3。

圖3 馬氏鏈的平穩(wěn)分布圖

從圖3發(fā)現(xiàn)，從第20周開(kāi)始，三種情況下分布都開(kāi)始收斂。最為奇特的是，不論初始狀態(tài)是三種行情中的哪一種，最終都會(huì)收斂到平穩(wěn)概率分布π＝（0．625，0．312，0．063）。也就是說(shuō)，收斂行為與初始狀態(tài)無(wú)關(guān)，收斂行為主要是由概率轉(zhuǎn)移矩陣P決定的。

所有的MCMC方法均以馬氏鏈定理為理論基礎(chǔ)。根據(jù)馬氏鏈定理，滿(mǎn)足一定條件的馬氏鏈的轉(zhuǎn)移序列會(huì)形成一個(gè)具有平穩(wěn)分布的樣本。一個(gè)自然的想法是，如果想獲得某個(gè)分布（稱(chēng)為目標(biāo)分布）的樣本，是否可以設(shè)計(jì)一個(gè)馬氏鏈（其平穩(wěn)分布為目標(biāo)分布），然后采集這個(gè)馬氏鏈的轉(zhuǎn)移序列來(lái)作為目標(biāo)分布的隨機(jī)樣本呢？答案是肯定的。

至此可以看到，用MCMC方法進(jìn)行隨機(jī)模擬的關(guān)鍵是設(shè)計(jì)一條馬爾科夫鏈，使其平穩(wěn)分布與目標(biāo)分布一致，而紛繁眾多的MCMC方法之間的區(qū)別就在于從不同的角度出發(fā)設(shè)計(jì)馬氏鏈。

（三）Metropol－Hastings算法：“千島湖植物考察”

為了更加生動(dòng)形象地說(shuō)明Metropol－Hastings算法，用“千島湖植物考察”為例進(jìn)行說(shuō)明①本例受John K．Kruschke（2011）Doing Bayesian Data Analysis：A Tutorial with R and BUGS書(shū)中"總統(tǒng)候選人巡島演講拉選票"章節(jié)啟發(fā)。。

某旅游風(fēng)景點(diǎn)“千島湖”由大大小小上千個(gè)島組成，千島湖風(fēng)景區(qū)植物種類(lèi)非常豐富。植物學(xué)家要去各小島進(jìn)行植物考察研究，考察的時(shí)間安排是依據(jù)各島上植物種類(lèi)數(shù)來(lái)確定的，植物種類(lèi)數(shù)多的島嶼考察時(shí)間長(zhǎng)，植物種類(lèi)數(shù)少的島嶼考察時(shí)間短。但在考察之前，植物學(xué)家并不清楚各個(gè)島的植物種類(lèi)數(shù)，甚至不知道一共有多少個(gè)島，所以不能事先確定留在每一個(gè)島上的天數(shù)，也無(wú)法得知在各島考察時(shí)間占總考察時(shí)間的比率——考察時(shí)間分布。那么，如何幫助植物學(xué)家安排在各島的行程路線(xiàn)和考察時(shí)間？為了用圖形進(jìn)行直觀(guān)介紹，將島嶼簡(jiǎn)化為9個(gè)，并假設(shè)島嶼一字排開(kāi)（如圖4中a），植物學(xué)家來(lái)到某一島嶼后，不僅可以了解該島上的植物種類(lèi)數(shù)，還可以了解相鄰島的植物種類(lèi)數(shù)。

精通MCMC的統(tǒng)計(jì)學(xué)者向植物學(xué)家建議用M－H算法的思想，通過(guò)“兩步隨機(jī)試探法”來(lái)決定考察路線(xiàn)和停留時(shí)間。植物學(xué)家到達(dá)某島嶼后，由于島嶼是一字排開(kāi)，所以面臨三種選擇（三個(gè)狀態(tài)）：一是繼續(xù)停留在島嶼上（保持原狀態(tài)）；二是到臨近左面島嶼（轉(zhuǎn)移到前狀態(tài)）；三是到臨近右面島嶼（轉(zhuǎn)移到后狀態(tài)）。這三種狀態(tài)可以轉(zhuǎn)化為兩種隨機(jī)選擇：一是去或留，二是左或右。但M－H算法突破常理，用“逆向思維”來(lái)確定選擇過(guò)程，先確定左或右，再?zèng)Q定去或留（究其原因可能是先基于左或右島的情況，再與當(dāng)前島嶼比較，最后決定是走還是留）?！皟刹诫S機(jī)試探法”的具體實(shí)施過(guò)程如下：

第一步，隨機(jī)確定方向—向左還是向右：植物學(xué)家在某島嶼上考察時(shí)（假設(shè)在第5個(gè)島嶼），就“可能”計(jì)劃前往鄰近島嶼（第4個(gè)島或第6個(gè)島），到底是向左還是向右呢？植物學(xué)家通過(guò)隨機(jī)擲一枚均勻硬幣來(lái)決定，如果硬幣為正面，就計(jì)劃向左去第4個(gè)島；如果硬幣為反面，則計(jì)劃向右去第6個(gè)島。

第二步，基于第一步的信息，部分隨機(jī)確定意愿—出發(fā)還是停留：假設(shè)第一步硬幣為正面，根據(jù)第一步規(guī)則，植物學(xué)家應(yīng)計(jì)劃“可能”向左去第4個(gè)島，但“可能”去，也“可能”不去。如何確定呢？統(tǒng)計(jì)學(xué)者給出的策略是，從當(dāng)前島（第5個(gè)島）的管理者那里可以得到目標(biāo)島（左面第4個(gè)島）的植物種類(lèi)數(shù)，比較目標(biāo)島和當(dāng)前島的植物種類(lèi)數(shù)量，再?zèng)Q定是否出發(fā)。

如果目標(biāo)島的植物種類(lèi)數(shù)比所處當(dāng)前島的植物種類(lèi)數(shù)多，那么植物學(xué)家就一定前往目標(biāo)島進(jìn)行考察研究；如果目標(biāo)島的植物種類(lèi)數(shù)比所處當(dāng)前島的少，那么植物學(xué)家又需要“另外一步隨機(jī)試探法”來(lái)確定意愿—依概率出發(fā)或是停留。

什么是依概率出發(fā)呢？舉例說(shuō)明：假設(shè)目標(biāo)島有150種植物，而所處當(dāng)前島有200種植物，則前往目標(biāo)島的概率就為0．75（150／200），繼續(xù)留在當(dāng)前島的概率為0．25（1－0．75）。那么植物學(xué)家到底是出發(fā)還是停留呢？似乎還沒(méi)有確切答案，此時(shí)的隨機(jī)判斷方法是：在地上畫(huà)一條1米長(zhǎng)的線(xiàn)段，然后隨機(jī)向該線(xiàn)段內(nèi)投一個(gè)石子，如果石子落在0和0．75米之間，則出發(fā)去目標(biāo)島考察研究；如果石子落在0．75和1之間，則繼續(xù)留在當(dāng)前島。

以上植物學(xué)家的隨機(jī)試探方法本質(zhì)上就是馬氏鏈的轉(zhuǎn)移矩陣，那么該植物學(xué)家設(shè)計(jì)的馬式鏈的平穩(wěn)分布（各個(gè)島嶼的考察研究的時(shí)間分布）是否與其目標(biāo)分布（即各個(gè)島嶼的植物種類(lèi)數(shù)分布）相同呢？用R軟件來(lái)模擬結(jié)果。

圖4 a目標(biāo)分布 b行動(dòng)軌跡 c頻數(shù)分布

圖4中a顯示出各個(gè)島的相對(duì)植物種類(lèi)數(shù)，這也是植物學(xué)家在各個(gè)島嶼逗留時(shí)間的目標(biāo)分布。圖4中b顯示出植物學(xué)家按照前述“兩步隨機(jī)試探法”進(jìn)行了5 000次轉(zhuǎn)移，轉(zhuǎn)移軌跡如圖4b所示。注意圖4中b的移動(dòng)軌跡：第一天，植物學(xué)家在最中間的島即第5個(gè)島上；第2天它仍然停留在第5個(gè)島嶼上；第3天轉(zhuǎn)移到了第4個(gè)島嶼上；第4天轉(zhuǎn)移到了第3個(gè)島嶼上；第5天又回到了第4個(gè)島嶼上……圖4中的c顯示在各個(gè)島嶼植物考察研究的頻數(shù)直方圖，可以看出各島嶼被考察的相對(duì)頻數(shù)，與圖4中a顯示出的各個(gè)島的相對(duì)植物種類(lèi)數(shù)相似，兩圖走勢(shì)是一樣的，這也說(shuō)明植物學(xué)家設(shè)計(jì)的馬式鏈的平穩(wěn)分布與目標(biāo)分布相同。

上述過(guò)程反映了M－H算法抽取樣本的基本原理?！扒u湖植物考察”的例子只是M－H算法的一種特殊情況：目標(biāo)分布是一維離散分布，建議分布也是離散的，分別以0．5的概率建議向左或向右移動(dòng)。M－H算法更一般的是研究高維連續(xù)型分布，建議分布也更復(fù)雜，但實(shí)質(zhì)與此例是相同的。以下給出M－H算法的一般步驟：

不斷循環(huán)以上過(guò)程，為什么就得到了目標(biāo)分布的隨機(jī)樣本了呢？尤其是從一個(gè)與π（x）毫不相關(guān)的建議分布q（x′｜x）產(chǎn)生隨機(jī)數(shù)，然后經(jīng)過(guò)一個(gè)簡(jiǎn)單的判斷過(guò)程，最終竟然得到了服從目標(biāo)分布的隨機(jī)數(shù)！事實(shí)上，建議分布q（x′｜x）和接受概率α（xt，x′）共同扮演了馬氏鏈的轉(zhuǎn)移核（類(lèi)似于離散馬式鏈的轉(zhuǎn)移矩陣）的作用，根據(jù)馬氏鏈定理，可證明其會(huì)收斂于目標(biāo)分布π（x）（見(jiàn)圖5）。

圖5 Metropolis－Hastings采樣軌跡圖

圖5所示為對(duì)一個(gè)二元正態(tài)分布進(jìn)行MCMC隨機(jī)模擬的路徑（實(shí)線(xiàn)），其中實(shí)線(xiàn)表示的是接受建議分布產(chǎn)生的轉(zhuǎn)移建議，細(xì)虛線(xiàn)表示建議分布產(chǎn)生的轉(zhuǎn)移建議被拒絕。由此可見(jiàn)，建議分布q（x′｜x）產(chǎn)生的轉(zhuǎn)移建議就像脫韁的“野馬”，四處亂竄，而接受概率就像一條隨機(jī)出現(xiàn)的“無(wú)形韁繩”，時(shí)時(shí)牽拽著這頭野馬，讓它在目標(biāo)分布所形成的“水草豐美”的區(qū)域內(nèi)漫步。

五、MCMC方法、R軟件與現(xiàn)代貝葉斯分析

貝葉斯定理簡(jiǎn)捷直觀(guān)，貝葉斯推斷清晰自然，由于應(yīng)用領(lǐng)域廣泛，受到了眾多學(xué)者的青睞。除了先驗(yàn)分布的選擇之外，貝葉斯分析在發(fā)展過(guò)程中遭遇最大的瓶頸就是后驗(yàn)分布的計(jì)算，特別是面對(duì)復(fù)雜高維的問(wèn)題時(shí)，難以用解析的方法求得明確的后驗(yàn)分布。

計(jì)算技術(shù)的進(jìn)步和MCMC方法的發(fā)展，重新喚醒了人們對(duì)貝葉斯分析這個(gè)古典統(tǒng)計(jì)思想的認(rèn)識(shí)。Metropolis等學(xué)者另辟蹊徑，沒(méi)有沿用先通過(guò)復(fù)雜的數(shù)學(xué)分析方法積分得到后驗(yàn)分布，再求其期望值、中位數(shù)和置信區(qū)間這些統(tǒng)計(jì)特征值的傳統(tǒng)方法，而是充分利用現(xiàn)代計(jì)算技術(shù)，基于馬爾可夫理論，使用蒙特卡洛模擬方法，回避后驗(yàn)分布表達(dá)式復(fù)雜的計(jì)算，創(chuàng)造性地使用MCMC方法，直接對(duì)后驗(yàn)分布的獨(dú)立隨機(jī)樣本進(jìn)行模擬，再通過(guò)分析模擬樣本來(lái)獲得相關(guān)統(tǒng)計(jì)特征值信息。MCMC方法的發(fā)展，解開(kāi)了貝葉斯分析計(jì)算的“緊箍咒”，現(xiàn)代貝葉斯分析重獲春天，再次生機(jī)勃發(fā)。

在近年頗為流行的統(tǒng)計(jì)分析軟件R中，基于MCMC采樣的現(xiàn)代貝葉斯分析越來(lái)越方便，許多可以實(shí)現(xiàn)MCMC的R程序包如雨后春筍般出現(xiàn)。在R主站上搜索關(guān)于MCMC的程序包，目前共有103個(gè)，而相對(duì)于傳統(tǒng)的貝葉斯計(jì)算軟件包（如BUGS、WinBUGS、OpenBUGS等），這些程序包功能更加強(qiáng)大與靈活。R中MCMC相關(guān)程序包大致可以分為5類(lèi)：

1．用于一般統(tǒng)計(jì)模型擬合。比較核心的包有：mcmc包，可以自定義后驗(yàn)分布函數(shù)，并從中進(jìn)行M－H采樣；MCMCpack包，包含了許多經(jīng)典統(tǒng)計(jì)模型的貝葉斯分析工具；bayesSurv包，提供了進(jìn)行貝葉斯生存回歸分析的函數(shù)；DPpackage包，所包含的函數(shù)可以完成貝葉斯非參數(shù)、半?yún)?shù)分析；bayesm包，可以用于微觀(guān)經(jīng)濟(jì)研究中的各種貝葉斯分析，如線(xiàn)性回歸模型、多項(xiàng)logit模型、多項(xiàng)probit模型、密度估計(jì)等。

2．用于特殊統(tǒng)計(jì)模型的貝葉斯分析。例如MCMCglmm包，專(zhuān)門(mén)求解貝葉斯廣義線(xiàn)性混合模型；lmm包，擬合貝葉斯線(xiàn)性混合模型；BayesTree包，完成貝葉斯回歸樹(shù)分析。

3．將R與其他采樣器連接起來(lái)的程序包。R的優(yōu)勢(shì)在于方便靈活的統(tǒng)計(jì)分析，進(jìn)行超大規(guī)模的模擬運(yùn)算可能遜色于專(zhuān)業(yè)的采樣器（如WinBUGS，BUGS）。R提供了與各種專(zhuān)業(yè)采樣器連接起來(lái)的接口，并封裝在程序包中，可自由加載使用。在R中定義貝葉斯分析模型，然后利用接口傳遞給采樣器并在采樣器中進(jìn)行MCMC模擬，模擬結(jié)束后，再利用接口將結(jié)果傳遞到R中，在R中對(duì)模擬結(jié)果進(jìn)行分析。這類(lèi)包有：R2WinBUGS（與 WinBUGS的接 口）；BRUGS（與 OpenBUGS 的 接 口）；rjags、R2jags與runjags（與JAGS的接口）。

4．處理MCMC樣本的包。coda包是核心的處理MCMC樣本的包，可以用于收斂診斷與輸出分析，幾乎所有關(guān)于MCMC的包都會(huì)依賴(lài)于coda包，足見(jiàn)此包的重要性。

5．輔助學(xué)習(xí)貝葉斯分析課程的包。這是R最人性化的一面，不僅提供了分析問(wèn)題的工具，還是學(xué)習(xí)貝葉斯分析方法的良師益友。在有關(guān)貝葉斯分析的參考書(shū)目里，經(jīng)典的當(dāng)屬《Bayesian Data Analysis》（《貝葉斯數(shù)據(jù)分析》Gelman，2003），為了方便學(xué)習(xí)，該書(shū)對(duì)應(yīng)的包BayesDA可以從CRAN加載使用，該包涵蓋了書(shū)中所涉及的數(shù)據(jù)集和函數(shù)，對(duì)于貝葉斯愛(ài)好者來(lái)說(shuō)，這無(wú)異于一學(xué)習(xí)利器。類(lèi)似的包還有 LearnBayes（《Bayesian Computation with R》《基于R的貝葉斯計(jì)算書(shū)中的對(duì)應(yīng)包》Jim Albert，2009．），此包由淺入深，其中編寫(xiě)的rwmetrop（）和rgibbs（）兩個(gè)函數(shù)對(duì)于理解 M－H算法和Gibbs采樣很有幫助。

六、結(jié)論與啟示

MCMC方法產(chǎn)生于物理和原子彈工作的研究中，在計(jì)算物理學(xué)（如粒子輸運(yùn)計(jì)算、量子熱力學(xué)計(jì)算、空氣動(dòng)力學(xué)計(jì)算）等領(lǐng)域應(yīng)用廣泛，之后基于MCMC的現(xiàn)代貝葉斯分析在生物醫(yī)學(xué)、金融工程學(xué)和宏觀(guān)經(jīng)濟(jì)學(xué)等方面的應(yīng)用，也得到了蓬勃的發(fā)展［15－16］。

2013年被稱(chēng)為“大數(shù)據(jù)元年”，隨著信息技術(shù)的疾速發(fā)展，自然科學(xué)中基因測(cè)序高維數(shù)據(jù)和社會(huì)科學(xué)中互聯(lián)網(wǎng)海量數(shù)據(jù)對(duì)統(tǒng)計(jì)分析理論和方法提出了嚴(yán)峻的挑戰(zhàn)。統(tǒng)計(jì)學(xué)向何處去？“離數(shù)學(xué)越來(lái)越遠(yuǎn)，與信息科學(xué)越來(lái)越近”的“統(tǒng)計(jì)與信息”相結(jié)合的趨勢(shì)似乎不可阻擋；基于云計(jì)算的隨機(jī)模擬算法、貝葉斯統(tǒng)計(jì)與頻率統(tǒng)計(jì)相融合的大數(shù)據(jù)分析新范式也初見(jiàn)倪端。

MCMC方法的發(fā)展和現(xiàn)代貝葉斯分析的復(fù)興改變了我們解決問(wèn)題的思路，它意味著一種思維范式的轉(zhuǎn)換?，F(xiàn)代貝葉斯拓寬了分析的視野，即從樣本空間的“有限范圍”放眼到樣本外“無(wú)邊區(qū)域”；MCMC方法改變了計(jì)算的重點(diǎn)，即從“精確解析”到“算法模擬”。MCMC帶領(lǐng)我們進(jìn)入一個(gè)全新的統(tǒng)計(jì)世界，在那里“小樣本”變成了“大數(shù)據(jù)”，“準(zhǔn)確”也被高精度的“近似”所替代。

60年前，當(dāng)Ulam和約翰·馮·諾依曼發(fā)明蒙特卡洛方法時(shí)，他們無(wú)法想象目前在處理基因組學(xué)和氣候?qū)W的巨型數(shù)據(jù)集時(shí)使用MCMC方法；250年前，當(dāng)理查德·普萊斯宣讀貝葉斯論文時(shí)，他也無(wú)法預(yù)測(cè)18世紀(jì)的貝葉斯定理會(huì)成為21世紀(jì)Google計(jì)算的新力量（用Google去搜索“18世紀(jì)的貝葉斯定理成為Google計(jì)算的新力量”）；當(dāng)今的統(tǒng)計(jì)學(xué)界是否也會(huì)有一種方法將被60年后的學(xué)者采用？是否也將有一個(gè)定理會(huì)被250年后的同仁們所紀(jì)念？這個(gè)問(wèn)題似乎可以基于MCMC方法，用貝葉斯定理去分析。

［1］ 陳希孺．?dāng)?shù)理統(tǒng)計(jì)學(xué)簡(jiǎn)史［M］．長(zhǎng)沙：湖南教育出版社，2002，50－51．

［2］ 劉樂(lè)平，袁衛(wèi)．現(xiàn)代貝葉斯分析與現(xiàn)代統(tǒng)計(jì)推斷［J］．經(jīng)濟(jì)理論與經(jīng)濟(jì)管理，2004（6）．

［3］ Andrieu C，De Freitas N，Doucet A，et al．An Introduction to MCMC for Machine Learning［J］．Machine learning，2003（1）．

［4］ Robert C．Casella G．A Short History of Markov Chain Monte Carlo：Subjective Recollections from Incomplete Data［J］．Statistical Science，2011（1）．

［5］ Metropolis N，Ulam S．The Monte Carlo Method［J］．Journal of the American Statistical Association，1949（4）．

［6］ Metropolis N，Rosenbluth A W，Rosenbluth M N，et al．Equations of State Calculations by Fast Computing Machines［J］．Journal of Chemical Physics，1953（2）．

［7］ Hastings W K．Monte Carlo Sampling Methods Usng Markov Cins and their Applications［J］．Biometrika，1970（1）．

［8］ Peskun P H．Optimum Monte－Carlo Sampling Using Markov Chains［J］．Biometrika，1973（3）．

［9］ Geman S，Geman D．Stochastic Relaxation，Gibbs Distributions and the Bayesian Restoration of Images［J］．IEEE，Transactions on Pattern Analysis and Machine Intelligance，1984（6）．

［10］Tanner M，Wong W．The Calculation of Posterior Distributions by Data Augmentation［J］．Journal of the American Statistical Association，1987（2）．

［11］Gelfand A，Smith A．Sampling Based Approaches to Calculating Marginal Densities［J］．Journal of the American Statistical Association，1990（4）．

［12］Roberts G O，Rosenthal J S．Coupling and Ergodicity of Adaptive Markov Chain Monte Carlo Algorithms［J］．Journal of applied probability，2007（4）．

［13］Green P．Reversible Jump MCMC Computation and Bayesian Model Determination［J］．Biometrika，1995（4）．

［14］Hobert J P，Jones G L，Presnell B，et al．On the Applicability of Regenerative Simulation in Markov Chain Monte Carlo［J］．Biometrika，2002（4）．

［15］丁東洋，周麗莉．基于貝葉斯方法的信用評(píng)級(jí)模型構(gòu)建與違約概率估計(jì)［J］．統(tǒng)計(jì)與信息論壇，2010（9）．

［16］周麗莉，丁東洋．基于MCMC模擬的貝葉斯分層信用風(fēng)險(xiǎn)評(píng)估模型［J］．統(tǒng)計(jì)與信息論壇，2011（12）．