一、引言

馬爾可夫鏈蒙特卡洛方法(Markov Chain Monte Carlo)，簡稱MCMC。其產(chǎn)生于20世紀50年代早期，是在貝葉斯理論框架下，通過計算機進行模擬的蒙特卡洛方法(Monte Carlo)。該方法將馬爾可夫(Markov)過程引入到Monte Carlo模擬中，實現(xiàn)抽樣分布隨模擬的進行而改變的動態(tài)模擬，彌補了傳統(tǒng)的蒙特卡羅積分只能靜態(tài)模擬的缺陷。

Metropolis等人在1953年首次提出了基于馬爾可夫鏈的蒙特卡羅方法，即 Metropolis 算法，并在最早的計算機上編程實現(xiàn)。Metropolis算法是首個普適的采樣方法，并啟發(fā)了一系列 MCMC 方法，所以人們把它視為隨機模擬技術騰飛的起點。 Metropolis算法這篇論文[1]被收錄在《統(tǒng)計學中的重大突破》中，《Computing in Science & Engineering》嘗試列出了對20世紀科學與工程的發(fā)展和實踐影響最大的十種算法：

• Metropolis Algorithm for Monte Carlo

• Simplex Method for Linear Programming

• Krylov Subspace Iteration Methods

• The Decompositional Approach to Matrix Computations

• The Fortran Optimizing Compiler

• QR Algorithm for Computing Eigenvalues

• Quicksort Algorithm for Sorting

• Fast Fourier Transform

• Integer Relation Detection

• Fast Multipole Method

Metropolis Algorithm for Monte Carlo 被列為十大算法之首。 用于蒙特卡洛的Metropolis算法定義了一個收斂的馬爾可夫鏈，其極限就是所需的概率分布。Metropolis算法及其推廣算法已被稱為蒙特卡洛馬爾可夫鏈技術(MCMC)，因為這些算法模擬了一個馬爾可夫鏈，從極限分布中獲取抽樣。

搞自然科學研究的人，MCMC應該是基本裝備。所謂平生不識MCMC, 便稱英雄也枉然。

二、MCMC的兩大貼身護衛(wèi)

MCMC由兩個MC組成，即馬爾可夫鏈(Markov Chain ，也簡稱MC)和蒙特卡洛法（Monte Carlo Simulation，簡稱MC）。要弄懂MCMC的原理我們首先得搞清楚馬爾科夫鏈的原理和蒙特卡羅法。

1、馬爾可夫鏈及其平穩(wěn)分布

1）定義

馬爾可夫鏈的命名來自俄國數(shù)學家安德雷.馬爾可夫以紀念其首次提出馬爾可夫鏈和對其收斂性質(zhì)所做的研究。

馬爾可夫鏈是一組具有馬爾可夫性質(zhì)的離散隨機變量的集合。所謂馬爾可夫性質(zhì)是指某一時刻狀態(tài)轉(zhuǎn)移的概率只依賴于它的前一個狀態(tài)。

具體地，馬氏鏈的數(shù)學定義:

對概率空間內(nèi)的離散隨機變量集合

若隨機變量的取值都在可數(shù)集S內(nèi)：

且隨機變量的條件概率滿足如下關系

則隨機變量集合X被稱為馬爾可夫鏈，可數(shù)集S被稱為狀態(tài)空間，馬爾可夫鏈在狀態(tài)空間內(nèi)的取值稱為狀態(tài)。

馬氏鏈某一時刻狀態(tài)轉(zhuǎn)移的概率只依賴于它的前一個狀態(tài)，這樣定義，可以大大簡化模型的復雜度，因此馬氏鏈在很多時間序列模型中得到廣泛的應用，比如循環(huán)神經(jīng)網(wǎng)絡RNN，隱式馬爾科夫模型HMM等，當然MCMC也需要它。

馬爾可夫過程通常分為3類：

• 時間、狀態(tài)都是離散的，稱為馬爾可夫鏈。

• 時間連續(xù)、狀態(tài)離散的，稱為連續(xù)時間的馬爾可夫鏈。

• 時間連續(xù)、狀態(tài)連續(xù)的，稱為馬爾可夫過程。

可見馬氏鏈是時間、狀態(tài)都是離散的馬爾可夫過程。

2）馬氏鏈轉(zhuǎn)移概率矩陣和收斂現(xiàn)象

我們先來看馬氏鏈的一個具體的例子。社會學家經(jīng)常把人按其經(jīng)濟狀況分成3類：下層(lower-class)、中層(middle-class)、上層(upper-class)，我們用1、2、3分別代表這三個階層。社會學家們發(fā)現(xiàn)決定一個人的收入階層的最重要的因素就是其父母的收入階層。如果一個人的收入屬于下層類別，那么他的孩子屬于下層收入的概率是0.65，屬于中層收入的概率是0.28，屬于上層收入的概率是0.07。事實上，從父代到子代，收入階層的變化的轉(zhuǎn)移概率如圖2-1和圖2-2

圖2-1

圖2-2

設狀態(tài)空間S={1,2,3} (1、2、3分別代表人的收入階層：下層、中層、上層), 使用矩陣的表示方式，轉(zhuǎn)移概率矩陣記為:

假設當前這一代人處在下層、中層、上層的人的比例是概率分布向量:

那么他們的子女的分布比例將是

他們的孫子代的分布比例將是

第n代子孫的收入分布比例將是

假設初始概率分布為π0=[0.21,0.68,0.11]，則我們可以計算前n代人的分布狀況如下

我們發(fā)現(xiàn)從第7代人開始，這個分布就穩(wěn)定不變了，這個是偶然的嗎？我們換一個初始概率分布π0=[0.75,0.15,0.1]試試看，繼續(xù)計算前n代人的分布狀況如下:

我們發(fā)現(xiàn)，到第9代人的時候, 分布又收斂了。最為奇特的是，兩次給定不同的初始概率分布，最終都收斂到概率分布 π=[0.286,0.489,0.225]，也就是說收斂的行為和初始概率分布 π0無關。這說明這個收斂行為主要是由概率轉(zhuǎn)移矩陣P決定的。我們計算一下 P的n階矩陣為

我們發(fā)現(xiàn)，當 n 足夠大的時候，這個P的n階矩陣的每一行都是穩(wěn)定地收斂到π=[0.286,0.489,0.225]這個概率分布。這個收斂現(xiàn)象并非是我們這個馬氏鏈獨有的，而是絕大多數(shù)馬氏鏈的共同行為。這個性質(zhì)同樣不光是離散狀態(tài)，連續(xù)狀態(tài)時也成立。

3）馬氏鏈定理

馬氏鏈定理： 如果一個非周期馬氏鏈具有轉(zhuǎn)移概率矩陣P, 狀態(tài)空間S, i,j∈S且它的任何兩個狀態(tài)是連通的，那么

這就是馬氏鏈的收斂定理。

由該定理可以得出如下結(jié)論：

<1>

<2>

<3> π 是方程 πP=π的唯一非負解。其中,

π稱為馬氏鏈的平穩(wěn)分布。

請注意兩種表示方式的不同：

• 表示狀態(tài)空間S中的各個狀態(tài)1,2,j...的概率：

• 表示不同時刻的概率分布

這個馬氏鏈的收斂定理非常重要，所有的 MCMC方法都是以這個定理作為理論基礎的。定理的證明相對復雜，一般的隨機過程課本中也不給證明，所以我們就不用糾結(jié)它的證明了，直接用這個定理的結(jié)論就好了。我們對這個定理的內(nèi)容做一些解釋說明：

• 非周期的馬爾科夫鏈：這個主要是指馬爾科夫鏈的狀態(tài)轉(zhuǎn)化不是循環(huán)的，如果是循環(huán)的則永遠不會收斂。幸運的是我們遇到的馬爾科夫鏈一般都是非周期性的；

• 任何兩個狀態(tài)是連通的：這個指的是從任意一個狀態(tài)可以通過有限步到達其他的任意一個狀態(tài)，不會出現(xiàn)條件概率一直為0導致不可達的情況；

• 馬爾科夫鏈的狀態(tài)數(shù)可以是有限的，也可以是無限的。因此可以用于連續(xù)概率分布和離散概率分布；

• π通常稱為馬爾科夫鏈的平穩(wěn)分布。

我們僅證明定理的第二個結(jié)論

證明：

上式兩邊取極限就得到

證明完畢。

4) 馬氏鏈的平穩(wěn)分布與細致平衡條件

給定一個馬爾可夫鏈，若在其狀態(tài)空間S存在概率分布

則π是該馬爾可夫鏈的平穩(wěn)分布。

細致平衡條件（Detailed Balance Condition）：給定一個馬爾科夫鏈，概率分布π和概率轉(zhuǎn)移矩陣P，?i，j∈S, (S是狀態(tài)空間)，如果下面等式成立:

則此馬爾科夫鏈具有一個平穩(wěn)分布（Stationary Distribution)。

需要注意，這個細致平衡條件僅是一個充分條件，而不是必要條件，也就是說存在具有平穩(wěn)分布的馬爾科夫鏈不一定滿足此細致平衡條件。

證明如下：

由條件推出

得

所以此馬爾科夫鏈具有一個平穩(wěn)分布。

從初始概率分布π0(x)出發(fā)，我們在馬氏鏈上做狀態(tài)轉(zhuǎn)移，記Xi的概率分布為πi(x)，則有

由馬氏鏈收斂的定理, 概率分布πi(x)將收斂到平穩(wěn)分布 π(x)。假設到第n步的時候馬氏鏈收斂，則有

所以

也就是說無論我們的起始狀態(tài)是何種概率分布，在狀態(tài)轉(zhuǎn)移矩陣 P 的作用下，經(jīng)過足夠大的 n 步轉(zhuǎn)移之后，最后都收斂到一個平穩(wěn)分布，且該分布是唯一不變的。即平穩(wěn)分布僅由狀態(tài)轉(zhuǎn)移矩陣P決定。

設馬爾可夫鏈的狀態(tài)空間S={ω1，ω2，..., ωi，...ωn}（可以理解為樣本空間在隨機過程中的表達），如果我們從一個具體的初始狀態(tài) x0∈S（x0可能是ω2，也可能是ωi）開始,沿著馬氏鏈按照概率轉(zhuǎn)移矩陣做跳轉(zhuǎn)，那么我們可以得到一個轉(zhuǎn)移狀態(tài)序列 ，馬氏鏈終于露出了它的廬山真面：

盡管馬爾可夫鏈最終會收斂到所需的分布，但是初始樣本可能會遵循非常不同的分布，尤其是在起點位于低密度區(qū)域的情況下。結(jié)果，通常需要一個老化期，所以我們會舍棄前面的老化樣本，得到我們需要的平穩(wěn)樣本集

這就是平穩(wěn)分布 π的狀態(tài)樣本序列。最后平穩(wěn)分布的樣本集可能達到幾萬、幾十萬甚至更大的規(guī)模。統(tǒng)計各個狀態(tài)出現(xiàn)的頻次，歸一化后畫出統(tǒng)計直方圖，就可以模擬目標分布。

5）基于馬爾可夫鏈的采樣過程

1> 輸入馬爾可夫鏈狀態(tài)轉(zhuǎn)移矩陣P, 狀態(tài)空間S，設定狀態(tài)轉(zhuǎn)移收斂次數(shù)為n1, 迭代次數(shù)為n2;

2> 從任意簡單概率分布采樣得到初始狀態(tài)值x0∈S;

3> for t=0 to n1+n2-1: 從條件概率分布p(x|x(t))中采樣，得到樣本x = x(t+1). 去掉前面的老化樣本，最終的樣本集

即為我們需要的平穩(wěn)分布對應的樣本集。

這里大家會問：

• 如何進行簡單概率分布采樣和條件概率分布采樣？

簡單概率分布可選擇正態(tài)分布或均值分布，正態(tài)分布根據(jù)Box–Muller變換法或接受-拒絕采樣法采樣，均值分布可直接計算機采樣；條件概率分布可由輸入轉(zhuǎn)移矩陣P得到概率數(shù)據(jù)，然后根據(jù)離散隨機變量逆變換法采樣。這些采樣方法后面有詳細介紹。

• 怎么證明得到的樣本集就是平穩(wěn)分布對應的樣本集？

兩種采樣方法來理解：

第1種方法，從簡單概率分布中采樣 10000個的初始狀態(tài)樣本，用最簡單粗暴的方法對每一個初始樣本都執(zhí)行一次n(n>n1)步轉(zhuǎn)移，狀態(tài)之間的轉(zhuǎn)移概率從該馬爾科夫鏈的狀態(tài)轉(zhuǎn)移矩陣獲取，最后依如下概率落入狀態(tài)空間中的一個狀態(tài)。

當初始狀態(tài)樣本足夠多時，由大數(shù)定律，最終落入到平穩(wěn)馬氏鏈的各個狀態(tài)的數(shù)量所占的比例逼近各個狀態(tài)的平穩(wěn)分布概率。但是這樣一來計算量相當大，每次只利用了采樣過程中的最后一個狀態(tài)，且要運行 10000次進入平穩(wěn)狀態(tài)的轉(zhuǎn)移過程，數(shù)據(jù)的利用率很低，浪費也很大。

第2種方法，只需要對一個初始狀態(tài)進行狀態(tài)轉(zhuǎn)移，就可以實現(xiàn)整個采樣過程。我們把它進入穩(wěn)態(tài)的那一刻記作t0，當 t0 時刻的狀態(tài)向下一個時刻t1 轉(zhuǎn)移時，它依照穩(wěn)態(tài)分布概率

（由n-步轉(zhuǎn)移概率的性質(zhì)Chapman–Kolmogorov等式可知，n-步轉(zhuǎn)移矩陣是其之前所有轉(zhuǎn)移矩陣的連續(xù)矩陣乘法，證明見附錄1）

進入到狀態(tài) 1、狀態(tài)2 或者狀態(tài)j當中的任意一個。以此類推，對于后面的所有轉(zhuǎn)移時刻：t2、t3......tn，當前狀態(tài)都能依照穩(wěn)態(tài)分布中的概率進入到狀態(tài)空間的任何其它狀態(tài)。進行N次轉(zhuǎn)移后得到N個采樣結(jié)果，N個采樣結(jié)果中各個狀態(tài)的數(shù)量所占的比例和平穩(wěn)分布是一致的。這種方法巧妙地用時間軸的采樣分布等價替代了空間軸的采樣分布，省時省力。

圖2-3示意了進入平穩(wěn)期的采樣全過程（假設總共3個狀態(tài)）：

圖2-3

6）馬爾可夫鏈的性質(zhì)

馬爾可夫鏈的4個性質(zhì)：不可約性、常返性、周期性和遍歷性。

• 不可約是指當馬爾可夫鏈為任何狀態(tài)時，其訪問其余所有狀態(tài)的概率都是正的，即隨機變量可以在任意狀態(tài)間轉(zhuǎn)移。

• 常返性指給定一個狀態(tài)則經(jīng)過有限的時間它一定會重新回到這個狀態(tài)。

• 非周期指馬爾可夫鏈的狀態(tài)的遷移不會陷入小的閉環(huán)。

• 遍歷性指馬爾可夫鏈的的所有狀態(tài)互通且均為非周期的正常返狀態(tài)。

• 平穩(wěn)馬爾可夫鏈也被稱為齊次馬爾可夫鏈。即轉(zhuǎn)移概率不隨時間而變。

• 滿足細致平衡條件的馬氏鏈具有一個平穩(wěn)分布。

• 由極限定理可知，遍歷鏈是平穩(wěn)馬爾可夫鏈。

• 有限狀態(tài)馬爾可夫鏈轉(zhuǎn)移概率的極限分布一定是平穩(wěn)分布。

• 若初始概率是平穩(wěn)分布，則任意時刻的絕對概率分布等于初始分布，也即為平穩(wěn)分布。

7) 馬爾可夫鏈的作用

馬爾可夫鏈在處理隨機動力學時對問題建模的強大作用。

馬爾可夫鏈預測理論在教育學、經(jīng)濟學、金融投資、生物學、農(nóng)作物栽培、地質(zhì)災變，...特別是水資源科學中都得到了極為廣泛地應用。

例如，馬爾可夫鏈在排隊理論、統(tǒng)計數(shù)據(jù)、生物種群進化建模、隱馬爾可夫模型信息理論、語音識別、音字轉(zhuǎn)換、詞性標注、蒙特卡洛采樣等都有應用。

2、蒙特卡洛法

1）定義

蒙特卡洛(Monte Carlo，MC)方法, 又稱隨機抽樣法，統(tǒng)計試驗法或隨機模擬法。

蒙特卡羅法是由馮·諾伊曼等人在20世紀40年代為研制核武器的需要而首先提出并以摩納哥的賭城Monte?Carlo來命名的一種計算方法。

蒙特卡洛法是數(shù)學建模十大算法之一，是思想和技巧的藝術品。它的魔力在于，對于那些規(guī)模極大的問題，求解難度隨著問題的維數(shù)（自變量個數(shù)）的增加呈指數(shù)級別增長、出現(xiàn)所謂的“維數(shù)的災難”的問題，傳統(tǒng)方法無能為力，而蒙特卡洛方法卻可以獨辟蹊徑，基于隨機仿真的過程給出近似的結(jié)果。

2）主要步驟

蒙特卡洛法是一類通過隨機抽樣過程來求取最優(yōu)解的方法的總稱，如果建立蒙特卡洛模型的過程沒有出錯，那么抽樣次數(shù)越多，得到的答案越精確。蒙特卡洛法可以歸納為如下三個步驟：

<a> 構(gòu)造或描述概率過程，將欲求解問題轉(zhuǎn)化為概率過程。

對于本身就具有隨機性質(zhì)的問題，主要是正確描述和模擬這個概率過程，對于本來不是隨機性質(zhì)的確定性問題，比如計算定積分，就必須事先構(gòu)造一個人為的概率過程，它的某些參量正好是所要求問題的解。即要將不具有隨機性質(zhì)的問題轉(zhuǎn)化為隨機性質(zhì)的問題。

<b> 實現(xiàn)從已知概率分布抽樣。

構(gòu)造了概率模型以后， 按照這個概率分布抽取隨機變量。這一般可以直接由軟件工具包（見附錄2，3，4）調(diào)用，或抽取均勻分布的隨機數(shù)構(gòu)造。這成為實現(xiàn)蒙特卡洛方法模擬實驗的基本手段，也是蒙特卡洛方法被稱為隨機抽樣的原因。

<c> 通過樣本計算各類統(tǒng)計量，此類統(tǒng)計量就是欲求問題的解。

一般說來，構(gòu)造了概率模型并能從中抽樣后（即實現(xiàn)模擬實驗后），我們就要確定一個隨機變量，作為所要求的問題的解，我們稱它為無偏估計。建立各種估計量，相當于對模擬實驗的結(jié)果進行考察和登記，從中得到問題的解。

3）組成部分

• 概率密度函數(shù)(pdf)：必須給出描述一個物理系統(tǒng)的一組概率密度函數(shù);

• 隨機數(shù)產(chǎn)生器：能夠產(chǎn)生在區(qū)間[a,b]上滿足概率分布的隨機數(shù);

• 抽樣規(guī)則：隨機抽取服從給定的pdf的隨機變量;

• 模擬結(jié)果記錄：記錄一些感興趣的量的模擬結(jié)果;

• 誤差估計：必須確定統(tǒng)計誤差（或方差）隨模擬次數(shù)以及其它一些量的變化;

• 減少方差的技術：利用該技術可減少模擬過程中計算的次數(shù);

• 并行和矢量化：可以在先進的并行計算機上運行的有效算法。

4）優(yōu)點

• 對于具有統(tǒng)計性質(zhì)問題可以直接進行解決；

• 對于連續(xù)性的問題不必進行離散化處理；

• 能夠比較逼真地描述具有隨機性質(zhì)的事物的特點及物理實驗過程；

• 受幾何條件限制??；

• 方法的誤差與問題的維數(shù)無關；

• 收斂速度與問題的維數(shù)無關；

• 具有同時計算多個方案與多個未知量的能力；

• 誤差容易確定；

• 程序結(jié)構(gòu)簡單，易于實現(xiàn)。

5）應用

蒙特卡洛方法應用非常廣泛，其領域包括統(tǒng)計物理、天體物理、物理化學、數(shù)學、計算生物、統(tǒng)計學、經(jīng)濟學、金融證券、社會科學等等。

蒙特卡洛方法在解決物理和數(shù)學問題時時常被使用。當十分困難或不可能得到解析表達式，或不可能應用確定性算法時，該方法顯得尤其重要和最為有用。

蒙特卡洛方法主要用于兩個不同類別的問題：

• 確定性的數(shù)學問題

無理數(shù)、不規(guī)則圖形的面積、多重積分、求逆矩陣、解線性代數(shù)方程組、解積分方程、解某些偏微分方程邊值問題和計算代數(shù)方程組、計算微分算子的特征值等等。

• 隨機性問題

仿真、概率分布的生成（MCMC算法）。

應用之一 求無理數(shù)π

解題思路：

構(gòu)造一個單位正方形和一個單位圓的1/4，往整個區(qū)域內(nèi)隨機投入點，根據(jù)點到原點的距離判斷點是落在1/4的圓內(nèi)還是在圓外，從而根據(jù)落在兩個不同區(qū)域的點的數(shù)目，求出兩個區(qū)域的比值。如此一來，就可以求出1/4單位圓的面積，從而求出圓周率π。

解題步驟：

第1步，建模，構(gòu)造概率過程，選擇區(qū)域內(nèi)隨機投入點作為隨機變量；

第2步，選擇二維均勻分布Y = sqrt(h^2+v^2), h、v~U[0,1], 作為已知概率分布抽樣；

第3步，根據(jù)落入正方形樣本數(shù)和1/4圓樣本數(shù)，求出結(jié)果。

演示代碼：

應用之二 求定積分

例子1

解題思路如下：

假定隨機變量具有密度函數(shù)

根據(jù)數(shù)學期望公式得

構(gòu)造統(tǒng)計數(shù)

根據(jù)大數(shù)定律

最終求出了定積分。

例子2 求定積分

解題思路：

將原式變形為

解：采用蒙特卡洛方法， MATLAB模擬，

N=500;

x=unifrnd(0,2,N,1);

mean(2*exp(3*x.^2))

第一條語句設定了停止條件，共做N次Monte Carlo 模擬；

第二條語句實現(xiàn)了在積分區(qū)間[0,2]均勻產(chǎn)生N個隨機數(shù)；

第三條語句求樣本均值，實現(xiàn)蒙特卡洛計算方法的面積逼近。

應用之三 仿真中子逸出

圖4-1是一個中子穿過用于中子屏蔽的鉛墻示意圖。鉛墻的高度遠大于左右厚度。設中子是垂直由左端進入鉛墻，在鉛墻中運行一個單位距離然后與一個鉛原子碰撞。碰撞后，任意改變方向，并繼續(xù)運行一個單位后與另一個鉛原子碰撞。這樣下去，如果中子在鉛墻里消耗掉所有的能量或者從左端逸出就被視為中子被鉛墻擋住，如果中子穿過鉛墻由右端逸出就視為中子逸出。如果鉛墻厚度為5個單位，中子運行7個單位后能量耗盡，求中子逸出的幾率。

解題思路：

這個問題并不復雜，但不容易找到一個解析表達式。而用模擬仿真的方法求解卻可以有滿意的結(jié)果。解題步驟如下：

第1步，建模，構(gòu)造概率過程，選擇中子碰撞后的角度作為隨機變量；

第2步，選擇角度的均勻分布作為已知概率分布抽樣；

第3步，根據(jù)假設的總中子樣本數(shù)和逸出的中子樣本數(shù)，求出結(jié)果。

解：具體做法如下：

我們關心的是一次碰撞后，中子在水平方向行進了多少，所以行進方向是正負θ的結(jié)果是一樣的，我們就只考慮θ是正的情形。由于中子運行的方向θ是隨機的，我們用計算機抽取在0到π間均衡分布的隨機數(shù)，模擬1000000個中子在鉛墻里行進的情形，看看這些中子與鉛原子碰撞7次后，有多少超過了鉛墻的右端。

實現(xiàn)偽代碼

n = 1000000;m = 0;flag = 1;for i=1:n    x=1;    for j=1:7        angle = pi * rand;        x = x + cos(angle);        if x<0            k=0;            flag=0;        end    end    if x>5 & flag==1        k=1;    else         k=0;    end    m = m + k;    flag = 1;end

運行結(jié)果：m/n = 1.3%

我們運行程序得出逸出鉛墻的中子的可能性約為1.3%，有了這個數(shù)字，我們可以報告安全部門，如果數(shù)字不能達到安全要求，我們則要加厚鉛墻。

3、馬爾可夫鏈與蒙特卡洛方法的靈魂附體

對于給定的概率分布p(x)，我們希望能有便捷的方式生成它對應的樣本。由于馬氏鏈能收斂到平穩(wěn)分布，于是一個很的漂亮想法是：如果我們能構(gòu)造一個轉(zhuǎn)移矩陣為P的馬氏鏈，使得該馬氏鏈的平穩(wěn)分布恰好是p(x), 那么我們從任何一個初始狀態(tài)x0出發(fā)沿著馬氏鏈轉(zhuǎn)移, 得到一個轉(zhuǎn)移狀態(tài)序列

如果馬氏鏈在第n步已經(jīng)收斂了，于是我們就得到了平穩(wěn)分布p(x)的樣本集

這個絕妙的想法在1953年被 Metropolis等人想到了，為了研究粒子系統(tǒng)的平穩(wěn)性質(zhì)， Metropolis 考慮了物理學中常見的波爾茲曼分布的采樣問題，首次提出了基于馬氏鏈的蒙特卡羅方法，即Metropolis算法，并在最早的計算機上編程實現(xiàn)。

三、MCMC有什么用？

1、貝葉斯推斷中很多問題都是MCMC方法的應用

從貝葉斯(Bayesian)的觀點來看，模型中的觀測變量和參數(shù)都是隨機變量，因此，樣本x與參數(shù)θ的聯(lián)合分布可以表示為：

根據(jù)貝葉斯理論，可以通過樣本x的信息對θ的分布的進行更新，后驗概率為：

其中π(θ)先驗概率；f(x|θ)似然函數(shù)（即樣本的概率分布）；分母是標準化常量。后驗概率 = (似然函數(shù)*先驗概率)/標準化常量 = 標準似然度*先驗概率。

貝葉斯學習中經(jīng)常需要進行三種積分運算：

• 歸范化（normalization)

• 邊緣 化（marginalization)

• 數(shù)學期望（expectation)

1> 后驗概率計算中需要歸范化計算：

2> 如果有隱變量z∈Z，后驗概率的計算需要邊緣化計算：

3> 如果有一個函數(shù)g(θ)，可以計算該函數(shù)的關于后驗概率分布的數(shù)學期望：

此積分值為樣本x的函數(shù)，因此可以對g(θ)進行推斷。

該期望估計

當維數(shù)很高時，直接進行這三類積分是非常困難的， MCMC方法就是為此目的而誕生，為這些計算提供了一個通用的有效解決方案。

2、MCMC的作用總結(jié)

1）在貝葉斯推理中主要用于從后驗分布的“非標準化部分”中直接生成樣本集，避免復雜計算。

對很多貝葉斯推斷問題來說，有時候后驗分布過于復雜，使得積分沒有顯示結(jié)果，數(shù)值方法也很難應用；有時候需要計算多重積分（比如后驗分布是多元分布時）。這些都會帶來計算上的很大困難。這也是在很長的時期內(nèi)，貝葉斯統(tǒng)計得不到快速發(fā)展的一個原因。1990年代MCMC計算方法引入到貝葉斯統(tǒng)計學之后，一舉解決了這個計算的難題。MCMC方法的最新發(fā)展使計算大型分層模型成為可能，這些模型需要對數(shù)百到數(shù)千個未知參數(shù)進行積分?？梢哉f，近年來貝葉斯統(tǒng)計的蓬勃發(fā)展，特別是在各個學科的廣泛應用和MCMC方法的使用有著極其密切的關系。

2）從任意給定的概率分布或復雜分布中獲取隨機樣本集

MCMC算法是一個普適的采樣方法，它的出現(xiàn)帶來了隨機模擬技術的騰飛。在隨機變量生成技術中，MCMC是一種相當高級的方法，可以從一個非常困難的概率分布中獲得樣本。更出乎意料的是，可以用MCMC從一個未經(jīng)標準化的分布中獲得樣本，這來自于定義馬爾可夫鏈的特定方式，馬爾可夫鏈對這些歸一化因子并不敏感。

對高維（多個自變量）積分的數(shù)值近似計算是MCMC方法的一種重要的應用，很多統(tǒng)計問題，比如計算概率、矩（期望、方差等）都可以歸結(jié)為定積分的計算。

MCMC的精髓在于構(gòu)造合適的馬氏鏈。

MCMC是一種簡單有效的計算方法，在很多領域到廣泛的應用，如統(tǒng)計學、貝葉斯(Bayes)問題、計算機問題等。

四、MCMC生成的樣本集長啥樣？

MCMC背后的基本思想是構(gòu)造一個遍歷的馬爾可夫鏈，使得其平穩(wěn)的、收斂的極限分布成為人們所需要的目標概率分布。

設馬爾可夫鏈的狀態(tài)空間S={ω1，ω2，...ωn}（可以理解為樣本空間在隨機過程中的表達），從某個初始狀態(tài)x0出發(fā)，用構(gòu)造的馬爾可夫鏈隨機游走，產(chǎn)生狀態(tài)樣本序列:

舍棄前面的老化(不穩(wěn)定)的樣本集合

應用遍歷定理，確定正整數(shù)m和n, 得到平穩(wěn)分布的狀態(tài)樣本集合：

該樣本集可能達到幾萬、幾十萬、幾千萬甚至更大的規(guī)模。統(tǒng)計該狀態(tài)樣本集中對應的各個狀態(tài)ω1，ω2...出現(xiàn)的頻次，歸一化后畫成直方圖。如圖3-1所示，綠色部分就是狀態(tài)概率分布（即抽樣概率分布），紅色部分是目標概率分布。t趨近∞時，抽樣概率分布可以無限逼近目標概率分布。如圖4-1所示。

圖4-1

五、MCMC與隨機采樣法

平時我們接觸比較多的場景是，給定一堆樣本數(shù)據(jù)，求出這堆樣本的概率分布p(x)。而采樣剛好是個逆命題：給定一個概率分布p(x)，如何生成滿足條件的樣本？

1、為什么需要采樣樣本？

通過樣本可以計算各類統(tǒng)計量（如均值、方差、矩等），而此類統(tǒng)計量就是欲求問題的近似解。例如求積分，根據(jù)大數(shù)定律，相當于求某概率分布采樣樣本的均值。

2、采樣必須滿足的條件

1）足夠地隨機

蒙特卡羅模擬有一個危險的缺陷：如果必須輸入一個模式中的隨機數(shù)并不像設想的那樣是隨機數(shù)，而卻構(gòu)成一些微妙的非隨機模式，那么整個的模擬和預測結(jié)果都可能是錯的。

對于給定的概率分布來說，采樣必須是隨機的，這個非常重要。隨機采樣對于人來說，是不能也，非不為也。依靠人來采樣很難得到真正的隨機樣本，人為地通過概率分布函數(shù)計算得到一批樣本注定是徒勞的。

蒙特卡羅模擬要取決于可靠的無窮盡的隨機數(shù)目來源。產(chǎn)生隨機數(shù)要靠隨機數(shù)發(fā)生器，隨機數(shù)發(fā)生器主要包括真（物理性）隨機數(shù)發(fā)生器和偽隨機數(shù)發(fā)生器。

真隨機數(shù)生成器（True Random Number Generator, TRNG）是一種通過物理過程而不是計算機程序來生成隨機數(shù)字的設備。通常是基于一些能生成隨機的“噪聲”信號的微觀現(xiàn)象，如熱力學噪聲、光電效應和量子現(xiàn)象等，其物理過程在理論上是完全不可預測的，并且這種不可預測性要經(jīng)過實驗檢驗。真隨機數(shù)生成器通常由換能器、放大器和模擬數(shù)字轉(zhuǎn)換器組成。其中換能器用來將物理過程中的某些效果轉(zhuǎn)換為電信號，放大器及其電路用來將隨機擾動的振幅放大到宏觀級別，而模擬數(shù)字轉(zhuǎn)換器則用來將輸出變成數(shù)字，通常是二進制的零和一。通過重復采樣這些隨機的信號，一系列的隨機數(shù)得以生成。

硬件隨機數(shù)發(fā)生器通常每秒僅產(chǎn)生有限數(shù)量的隨機位。為了增加可用的輸出數(shù)據(jù)速率，通常將它們用于生成“種子”，以用于更快的密碼安全偽隨機數(shù)生成器，然后生成器以更高的數(shù)據(jù)速率生成偽隨機輸出序列。

偽隨機數(shù)發(fā)生器是依靠計算機按照一定算法模擬產(chǎn)生的，其結(jié)果是確定的，是可見的。通過這種方式產(chǎn)生的“隨機數(shù)”并不隨機，是偽隨機數(shù)。貝爾實驗室的里德博士告誡人們記住偉大的諾依曼的忠告:“任何人如果相信計算機能夠產(chǎn)生出真正的隨機的數(shù)序組都是瘋子?！?/p>

偽隨機數(shù)存在兩大問題

• 遞推公式和初始值確定后，整個隨機數(shù)序列便被唯一確定，不滿足隨機數(shù)相互獨立的要求。

• 隨機數(shù)序列會出現(xiàn)周期性的循環(huán)現(xiàn)象。

真隨機數(shù)的生成對技術的要求比較高，在實際應用中往往使用偽隨機數(shù)就夠了。在使用偽隨機數(shù)時，要盡量克服其存在的兩大問題。對于第一個問題，不能從本質(zhì)上加以改變，但只要遞推公式選的比較好，隨機數(shù)間的相互獨立性是可以近似滿足的；對于第二個問題，因為用蒙特卡羅方法解任何具體問題時，所使用的隨機數(shù)的個數(shù)總是有限的，只要所用隨機數(shù)的個數(shù)不超過偽隨機數(shù)序列出現(xiàn)循環(huán)現(xiàn)象時的長度就可以了。

2）足夠多的采樣樣本

我們通過采樣獲取足夠多的樣本， 然后統(tǒng)計樣本空間中的樣本出現(xiàn)的頻次，來模擬和逼近目標概率分布的。采樣樣本越多，越逼近目標概率分布。如圖5-1所示。

幾萬或幾十萬的樣本集是小目標；幾百萬或幾千萬的樣本集很常見；在AI和大數(shù)據(jù)時代，目標概率分布有時會涉及到幾百個、幾千個參數(shù)， 幾百億的樣本集都出現(xiàn)了。

圖5-1

3） 盡量覆蓋樣本空間

這與算法和采樣策略有關。

4） 能模擬給定的目標分布

這與算法有關，要從理論上證明。

3、常見的隨機采樣方法

1） U~[0，1均勻分布采樣法（Uniform Sampling）

例如：線性同余偽隨機數(shù)生成器

X是偽隨機序列；

m, m>0表示模量；

a, 0<a<m表示乘數(shù)；

c, 0≤c<m表示增量；

X0， 0≤X0<m表示初始值；

如果c=0, 稱為乘法同余發(fā)生器，如果c!=0, 稱為混合同余發(fā)生器；

Yn則是 [0, 1)內(nèi)服從均勻分布的隨機變量。

小結(jié)：

均勻分布隨機數(shù)可以直接由計算機模擬產(chǎn)生。

2） Box–Muller變換法

Box–Muller 變換是一種利用均勻分布快速產(chǎn)生符合標準正態(tài)分布偽隨機數(shù)對的一種方法。基本思想是先得到服從均勻分布的隨機數(shù)，再將服從均勻分布的隨機數(shù)轉(zhuǎn)變?yōu)榉恼龖B(tài)分布。

I 標準形式的 Box–Muller 變換

假設變變量 U1 和變量 U2 是（0，1]均勻分布的隨機數(shù)，

且 U1 和 U2 彼此獨立，令：

則 Z0 和 Z1 就是服從 N（0，1）的標準正態(tài)分布隨機數(shù)，并且 Z0 和 Z1 相互獨立。

II 極坐標形式的 Box–Muller 變換

假設變量 u 和變量 v 是在[-1，1]上的均勻分布隨機量，

u 和 v 相互獨立，令：

故而，隨機數(shù) z0 和 z1 計算后得出如下結(jié)果：

z0 和 z1 是服從分布 N（0，1）的隨機數(shù)，并且 z0 和 z1 相互 獨立。

3） 逆變換法（Inverse Transform Method）

逆變換法就是利用累積分布函數(shù)求逆，將所求概率分布的抽樣轉(zhuǎn)化為隨機均勻分布的抽樣。

I 連續(xù)隨機變量逆變換法

設所求連續(xù)分布采樣的連續(xù)隨機變量X的概率密度函數(shù)為f (x) 。

第1步，求概率累積分布函數(shù)F(x)

F(x)是單調(diào)遞增的連續(xù)函數(shù)，令

則U也是隨機變量，且U∈[0, 1]。

第2步，求反函數(shù)

第3步，隨機產(chǎn)生均勻分布U[0,1]，求出對應的X，就可得到所求概率分布的采樣。

II 離散隨機變量逆變換法

設所求離散分布采樣的離散隨機變量的概率密度為

第1步，生成均勻分布U[0,1]隨機數(shù)。

第2步，求出對應的X

也就是如下分段式：

X就是所求概率分布的采樣。

III 輪盤賭采樣

任意離散分布都可以畫在輪盤上，采樣只需要隨機地旋轉(zhuǎn)輪盤即可。如圖5-2所示。

該輪盤模擬了累積分布函數(shù)，隨機轉(zhuǎn)動輪盤相當于隨機均勻采樣，最后輪盤的停止位置就是采樣值。

圖5-2

小結(jié)：

優(yōu)點：隨機數(shù)可以直接由均勻分布導出。

缺點：很多復雜累積分布函數(shù)求逆困難。

4）接受-拒絕采樣法（Acceptance-Rejection Sampling）

實際隨機分布抽樣中會遇到如下問題：

• 概率密度函數(shù)p(x) 太復雜無法直接采樣；

• 逆變換中一些累積概率分布函數(shù)無法求逆。

接受-拒絕采樣法就派上了用場?；舅枷胧沁x擇一個容易抽樣的隨機分布函數(shù)q(x)和常數(shù)k, 使得k*q(x)≥p(x)（這里p(x)可以進行非歸一化，乘以常數(shù)Z，以下省略），然后按照一定的方法拒絕某些抽樣，達到接近p(x)分布的目的。

如圖5-3所示：

圖5-3

假設我們選取已知的高斯分布q(x)=Gassian(1.4, 1.2)，我們按照一定的方法拒絕某些樣本，達到接近p(x)分布的目的。

具體步驟如下：

第1步，選擇已知的分布函數(shù)q(x)，確定常量k，使得p(x)總在kq(x)的下方。

第2步，在x軸上，從q(x)分布抽樣取得a。但是a不一定留下，會有一定的幾率被拒絕。

第3步，在y軸上，從均勻分布(0,1)中抽樣得到u。如果u.kq(a)>p(a)，就是落到了灰色的區(qū)域中，拒絕；否則接受這次抽樣。

第4步，重復2、3步。

模擬結(jié)果如圖5-4所示：

圖5-4

小結(jié)：

在高維（多個自變量）分布的情況下，接受-拒絕采樣會出現(xiàn)如下一些問題：

1> 合適的提議分布q(x)比較難以找到，并且當kq(x)與p(x)相差太多時，會導致拒絕率很高，無用計算增加。

2> 很難確定一個合理的常數(shù) k 值, 并且k值往往會很大；

3> 所需的樣本量隨空間維數(shù)增加而指數(shù)增長。

5） 重要性采樣法（Importance Sampling）

重要性采樣是統(tǒng)計學中估計某一分布性質(zhì)時使用的一種方法，該方法從與原分布不同的新分布中采樣，而對原分布的性質(zhì)進行估計。重要性采樣常用于蒙特卡洛積分。

作用一 降方差

重要性采樣是蒙特卡洛方法中的一個重要策略。該方法不改變統(tǒng)計量，只改變概率分布。重要性采樣算法就是在有限的采樣次數(shù)內(nèi)，盡量讓采樣點覆蓋對積分貢獻很大的點或者以一種受控的方式改變仿真，以便增加稀少事件的數(shù)目。它是一種降方差的仿真方法。

作用二 將復雜的無法抽樣的分布轉(zhuǎn)化為簡單的容易抽樣的分布，間接求一些統(tǒng)計量f(x)

其中， q(x)為提議分布(proposal distribution)，p(x) / q(x)可以看做 importance weight。

小結(jié)：

優(yōu)點： 如果我們無法從分布p(x)中抽樣, 或者從中抽樣的成本很高, 那么我們還可以求E[f(x)]嗎？答案是可以的, 我們可以從一個簡單的分布q(x)中進行抽樣得到x, 求f(x)p(x)/q(x)的期望，進而可以求得函數(shù)積分的近似解。

缺點：在高維空間里找到一個這樣合適的q(x)非常難。

6） MCMC采樣法

（1）MCMC思想的來源

拒絕采樣法和重要性采樣法面臨兩個問題：

1> 如果提議分布q(x)與目標分布p(x)差距過大，模擬采樣效果就不太好；

2> 構(gòu)造一個與目標分布p(x)相似的提議分布q(x)很困難。

我們知道，對于任意的目標分布，我們可以設計一種接受-拒絕的概率，從而使得所有接受的點都是該分布的樣本。拒絕采樣法和重要性采樣法采用了提議函數(shù)q(x)和接受概率實現(xiàn)了這種思想，然而問題在于，由于其提議分布q(x)是固定的，這對q(x)的要求很高，試想下，如果q與真實p的差距過大，那么幾乎每個樣本都會被拒絕掉，效率很低。那如何設計更好的q呢？一個巧妙的辦法是設計一個自適應的提議分布 q(x’|x)，它以上一次接受的樣本作為條件，然后經(jīng)過轉(zhuǎn)移propose一個新的樣本，如果我們能夠約束這個轉(zhuǎn)移不會改變目標分布，那么我們就能快速的找到目標分布p(x)的樣本。如圖5-5所示，這也正是MCMC的思想。

圖5-5

另一方面，大多數(shù)拒絕采樣法和重要性采樣法都遭受“維數(shù)的詛咒”，其中拒絕的概率隨維數(shù)的增加呈指數(shù)增長。 MCMC方法都沒有這個“維數(shù)詛咒”的問題，當要采樣的分布維數(shù)很大時，它通常是唯一可用的解決方案。結(jié)果，MCMC方法成為了從分層貝葉斯模型和當今許多學科中使用的其它高維統(tǒng)計模型中生成樣本的首選方法。MCMC采樣是一種普適采樣方法，它可以解決其它采樣方法無法搞定的很多問題。

（2）MCMC中的提議分布q(.)

提議分布q(x, x’) = q(x’|x)通常要滿足以下三個條件：

1> 對于固定的x, q(.)是一個概率密度函數(shù);

2> 對于?x, x’∈S, q(x, x’)的值要能計算;

3> 對于固定的x, 能夠方便地從q(x,x’)中產(chǎn)生隨機數(shù)。

從理論上講，提議函數(shù)q(.)的選取是任意的，但在實際計算中，提議函數(shù)的選取對于算法效率的影響是相當大的。一般認為提議函數(shù)與目標分布越接近，模擬的效果越好。

q(.)函數(shù)的選擇一般有兩種：

1> 對稱: 如高斯分布，均勻分布，這種情況下稱為隨機游走;

2> 非對稱：如 log-normal，采樣值傾向于選擇高概率的點，因為我們不希望采樣點來回震蕩。

（3） MCMC采樣原理

MCMC采樣法的關鍵是構(gòu)建馬氏鏈，馬氏鏈的收斂性質(zhì)主要由狀態(tài)轉(zhuǎn)移矩陣P決定, 所以基于馬氏鏈做采樣的關鍵問題是如何構(gòu)造狀態(tài)轉(zhuǎn)移矩陣P, 使得馬氏鏈具有一個平穩(wěn)分布，且該平穩(wěn)分布π(x)恰好是我們要的分布p(x)。如何能做到這一點呢？

答案的關鍵在細致平衡條件

其中i，j∈S (S是狀態(tài)空間)，P是狀態(tài)轉(zhuǎn)移矩陣。

不幸的是，對于平穩(wěn)分布π(x)，即我們要的分布p(x)，要找到滿足細致平衡條件的狀態(tài)矩陣P很困難。沒關系，我們隨機找一個已知的馬爾科夫鏈狀態(tài)轉(zhuǎn)移矩陣Q，其中Q(i,j)=Q(j|i)是從狀態(tài)i轉(zhuǎn)移到狀態(tài)j的提議分布條件概率，i,j∈狀態(tài)空間S，它可能一開始并不滿足細致平衡條件，即：

我們可以對上式做一個改造，使細致平衡條件成立。方法是引入一個α(i,j) ,使上式可以取等號，即：

要使該等式成立，只要滿足下面兩式即可

這樣，我們就得到了平穩(wěn)分布π(x)對應的馬爾科夫鏈狀態(tài)轉(zhuǎn)移矩陣P

我們有一般稱α(i,j)為接受概率

其取值在(0,1)之間，可以理解為一個概率值。

也就是說，基于平穩(wěn)分布的狀態(tài)轉(zhuǎn)移矩陣P可以通過任意一個馬爾科夫鏈狀態(tài)轉(zhuǎn)移矩陣Q以一定的接受概率獲得。在接受-拒絕采樣法中，那里是以一個常見分布通過一定的接受-拒絕概率得到目標概率分布，這里是以一個常見的馬爾科夫鏈狀態(tài)轉(zhuǎn)移矩陣Q通過一定的接受-拒絕概率得到目標轉(zhuǎn)移矩陣P，兩者解決問題的思路是類似的。

由于α(i,j) =π(j)Q(j, i)<1, 通常是一個比較小的數(shù)，因此在一次采樣中很容易拒絕跳轉(zhuǎn)，從而導致采樣的時間成本、計算成本很高。

因為α(i,j)<1和α(j,i)<1, 我們可以將兩者同時放大，且不破環(huán)等式成立。

1> 當α(i,j) ≤ α(j,i)時，則α(i,j)按照α(j,i)放大到1（概率不可能大于1）的比例等比例放大:

2> 當α(i,j) > α(j,i)時，則α(j,i)按照α(i,j)放大到1的比例等比例放大:

所以將兩者同時放大有

此時α(i,j)就是Metropolis-Hasting算法的接受概率。其中π(x)=p(x), Q(i,j)=Q(j|i)是提議分布條件概率，α(i,j)在取到1的情況下能實現(xiàn)鏈的滿概率跳轉(zhuǎn)，否則以放大1/(π(i)Q(i,j))倍的概率進行跳轉(zhuǎn)。

當提議分布對稱時，即Q(j, i) = Q(i, j)， 接受概率

此時α(i,j)就是Metropolis算法的接受概率。其中π(x)=p(x)。

由任選的提議分布Q(i,j)和接受概率α(j,i), 可以構(gòu)造一個滿足細致平衡條件的馬氏鏈， 該馬氏鏈具有一個平穩(wěn)分布。構(gòu)造平穩(wěn)分布的馬氏鏈的過程就是我們實現(xiàn)MCMC采樣的過程。

（4）MCMC幾種基本算法

MCMC算法包括Metropolis算法、Metropolis-Hasting 采樣法、吉布斯采樣法等。

I、 Metropolis算法

已知目標概率分布p(x)，如何對其進行采用？

我們分別采用接受-拒絕算法和Metropolis算法進行采樣，并對兩者做了對比。

接受-拒絕算法的具體步驟如下：

1> 首先選擇一個容易抽樣的分布作為提議分布，記為q(x)，接著確定一個常數(shù)C，其中 C > 1， 使得在x的定義域上都有p(x) ≤ Cq(x)。

2> 生成服從提議概率密度函數(shù)q(x)的隨機數(shù)y。

3> 接著生成一個服從均勻分布U (0,1) 的隨機數(shù)u。

4> 計算接受概率

如果u < α(y)，則接受y; 否則，拒絕。

5> 不斷重復2~ 4就可以生成一組服從目標分布p(x)的隨機數(shù)序列。

Metropolis算法就是改進的接受-拒絕算法，具體步驟如下：

1> 首先構(gòu)造提議分布q(.|x(t)), 這里提議分布是對稱的，即q(y|x)=q(x|y)， 例如，我們可以選均值為u=x(t)的正態(tài)分布為提議分布。設置t=0和初始狀態(tài)x0。

2> 從提議概率分布q(.|x(t))中產(chǎn)生一個隨機數(shù)y作為下一個狀態(tài)。

3> 接著生成一個服從均勻分布U (0,1) 的隨機數(shù)u。

4> 計算接受概率

如果u ≤ α，則接受候選隨機數(shù)y, 賦值x(t+1)=y，y作為提議分布的新狀態(tài)，意味著此操作將鼓勵跳轉(zhuǎn)到目標分布的高密度區(qū)。

如果u > α，則拒絕候選隨機數(shù)，賦值x(t+1)=x(t)，表明f(x)在此y狀態(tài)為低密區(qū), 以概率的形式被拒絕。

5> 增量設置t = t+1。

6> 重復步驟2~步驟5，我們就可以得到一組僅依賴于上一個隨機數(shù)并且與前面隨機數(shù)無關的隨機序列

Metropolis算法與接受-拒絕算法不同之處在于：

1> Metropolis算法的提議分布是自適應的條件概率分布，產(chǎn)生的隨機數(shù)是由上一次的狀態(tài)值轉(zhuǎn)移得到；接受-拒絕算法的提議分布是固定的，沒有狀態(tài)轉(zhuǎn)移的概念。

2> 兩者的接受概率定義不同。

Metropolis算法中唯一的限制是提議分布必須是對稱的，滿足條件的這類分布有正態(tài)分布、柯西分布、t分布和均勻分布等。為了能夠使用非對稱分布，或者目標分布的支持集是不對稱的，如x∈(0,+∞), 我們將考慮使用Metropolis-Hastings算法。

例子已知： (貝葉斯推斷: 一個簡單的投資模型) 假設有5種股票被跟蹤記錄了250的交易日每天的表現(xiàn), 在每一個交易日, 收益最大的股票被標記出來。用Xi表示股票i在250個交易日中勝出的天數(shù), 則記錄到得的頻數(shù)(x1,x2,... , x5)為隨機變量(X1,X2,...,X5) 的觀測?；跉v史數(shù)據(jù), 假設這5種股票在任何給定的一個交易日能勝出的先驗機會比率為1 : (1 ? β) : (1 ? 2β) : 2β : β, 這里β∈(0, 0.5)是一個未知的參數(shù)。在有了當前這250個交易日的數(shù)據(jù)后, 使用Bayes方法對此比例進行更新。

求：根據(jù)觀測的數(shù)據(jù)求β目標分布隨機抽樣的樣本集，進而估計β？

解：根據(jù)前述，(X1,X2,...,X5)在給定β的條件下服從多項分布, 概率向量為

因此β服從后驗分布為（根據(jù)聯(lián)合密度的多項分布公式，把x1,x2,...,x5看作常數(shù)）

我們不能直接從此后驗分布中產(chǎn)生隨機數(shù)。一種估計β的方法是產(chǎn)生一個鏈, 使其平穩(wěn)分布為此后驗分布, 然后從產(chǎn)生的鏈中抽樣來估計β.

我們這里使用隨機游動Metropolis算法, 可以初始化β=0.2, 提議分布為均勻分布U[0,1], 收斂次數(shù)burn=1000, 迭代次數(shù)m=5000, 接受的概率為

因為后驗分布即為所求的目標分布f(x), 可得

觀測數(shù)據(jù)(x1, x2, x3, x4, x5)是5只股票服從概率向量p在250天中勝出的頻次。

R語言的演示代碼：

小結(jié)：

Metropolis算法相鄰樣本是相關的，無法正確反映分布。這意味著，如果我們想要一組獨立的樣本，最后還必須丟棄一些樣本?？梢酝ㄟ^增加跳躍寬度選取樣本來降低自相關。

II、 Metropolis-Hastings算法

Metropolis算法的應用受到了提議概率分布的對稱形式的限制，因此，Hastings把提議概率分布從對稱分布形式推廣到一般情況，形成了Metropolis-Hastings算法，簡稱M-H算法，也稱之為通用梅特羅波利斯算法。

<1>、M-H算法對于一維目標概率分布的采樣

M-H算法具體步驟如下：

1> 首先構(gòu)造提議概率分布q(.|x(t))。設置t=0和初始狀態(tài)x0。

2> 從提議分布q(.|x(t))中產(chǎn)生一個隨機數(shù)y作為下一個狀態(tài)。

3> 接著生成一個服從均勻分布U (0,1) 的隨機數(shù)u。

4> 計算接受概率

如果u ≤ α，則接受候選隨機數(shù)y, 賦值x(t+1)=y;

如果u > α，則拒絕候選隨機數(shù)，賦值x(t+1)=x(t)。

5> 增量設置t = t+1。

6> 重復步驟2~步驟5，我們就可以得到一組僅依賴于上一個隨機數(shù)并且與前面隨機數(shù)無關的隨機序列。

采樣演示如下：

第1次采樣，選擇初始狀態(tài)x0

圖5-6a

第2次采樣，接受候選樣本x1

第3次采樣，接受候選樣本x2

圖5-6c

第4次采樣，拒絕候選樣本x3

圖5-6d

第5次采樣，接受候選樣本x4

圖5-6e

<2>、M-H算法對于高維目標概率分布的采樣

滿條件分布定義

馬爾可夫鏈蒙特卡羅法的目標分布通常是多元聯(lián)合概率分布

如果條件概率分布

中所有 k個變量全部出現(xiàn)，其中

那么稱這種條件概率分布為滿條件分布（full conditional distribution)。

滿條件分布有以下性質(zhì)：

1> 對任意的x, x’和任意的I?K, 有

2> 對任意的x, x’和任意的I?K, 有

M-H算法中，可以利用上述性質(zhì)，簡化計算，提高計算效率。具體地，通過滿條件分布概率的比

代替計算聯(lián)合概率的比

而前者更容易計算。

單分量Metropolis- Hastings算法

在Metropolis-Hastings算法中，通常需要對多元變量分布進行抽樣，有時對多元變量分布的抽樣是困難的。

可以對多元變量的每一變量的條件分布依次分別進行抽樣，從而實現(xiàn)對整個多元變量的一次抽樣，這就是單分量Metropolis-Hastings (single-component Metropolis-Hastings）算法。

假設馬爾可夫鏈的狀態(tài)由k維隨機變量表示

馬爾可夫鏈在時刻i的狀態(tài)

為了生成容量為的樣本集合

單分量Metropolis-Hastings 算法的基本做法：假設有一個k維的隨機向量，現(xiàn)想要構(gòu)造一條k維向量n個樣本的馬爾可夫鏈序列，那么隨機初始化一個k維向量，然后固定這個向量其中的k-1個元素，抽取剩下的那個元素，生成轉(zhuǎn)移狀態(tài)的隨機數(shù)，這樣循環(huán)k次，就把整個向量更新了一遍，也就是生成了一個新的k維向量樣本，把這個整體重復n次就得到了一條馬爾科夫鏈。

單分量M-H算法由下面的k步迭代實現(xiàn)M-H算法的一次迭代

馬爾可夫鏈的轉(zhuǎn)移概率為

如圖5-7所示，單分量M-H算法的迭代過程。目標是對含有兩個變量的隨機變量x進行抽樣

如果變量x1或x2更新，那么在水平或垂直方向產(chǎn)生一 個移動，連續(xù)水平和垂直移動產(chǎn)生

一個新的樣本點。

圖5-7

例子 邏輯回歸應用

考慮54位老年人的智力測試成績(Wechesler Adult Intelligence Scale, WAIS, 0-20分)。 研究的興趣在于發(fā)現(xiàn)老年癡呆癥。

我們采用如下簡單的邏輯回歸模型:

則似然函數(shù)為

考慮β0, β1的先驗分布π為獨立的正態(tài)分布:

從而我們需要從此分布中產(chǎn)生隨機數(shù)。

我們考慮如下三種抽樣方法.

M-H算法: 獨立抽樣

提議分布取為

則算法如下：

R語言演示代碼

wais-read.table (.wais. txt. ,header=TRUE)y<-wais[,2]; x<-wais[,1]m <- 55000 #length of chainmu.beta<-c(0,0); sigma.beta<-c(100,100)prop.s<-c(0.1,0.1) #proposal distribution standard variancebeta <- matrix(nrow=m, ncol=2)acc.prob <- 0current.beta<-c(0,0)for (t in 1:m){    prop.beta<- rnorm(2, current.beta, prop.s )    cur.eta<-current.beta[1]+current.beta[2]*x    prop.eta<-prop.beta[1]+prop.beta[2]*x    loga <-(sum(y*prop.eta-log(1+exp(prop.eta)))    -sum(y*cur.eta-log(1+exp(cur.eta)))    +sum(dnorm(prop.beta, mu.beta,s.beta,log=TRUE))    -sum(dnorm(current.beta,mu.beta,s.beta,log=TRUE)))    u<-runif(1)    u<-log(u)    if( u < loga) {        current.beta<-prop.beta        acc.prob <- acc.prob+1    }    beta[t,]<-current.beta}acc.prob<-acc.prob/macc.prob# convergence diagnostics ploterg.mean<-function( x ){     # compute ergodic mean    n<-length(x)    result<-cumsum(x)/cumsum(rep(1,n))}burnin<-15000idx<-seq(1,m,50)idx2<-seq(burnin+1,m)par(mfrow=c(2,2))plot(idx,beta[idx,1],type='l',xlab='Iterations',ylab='Values of beta0')plot(idx,beta[idx,2],type='l',xlab='Iterations',ylab='Values of beta1')ergbeta0<-erg.mean(beta[,1])ergbeta02<-erg.mean(beta[idx2,1])ylims0<-range(c(ergbeta0,ergbeta02))ergbeta1<-erg.mean(beta[,2])ergbeta12<-erg.mean(beta[idx2,2])ylims1<-range(c(ergbeta1,ergbeta12))plot(idx , ergbeta0[idx], type=’l’, ylab=’Values of beta0’, xlab=’Iterations’,main=’(c) Ergodic Mean Plot of beta0’, ylim=ylims0)lines(idx2, ergbeta02[idx2-burnin], col=2, lty=2)plot(idx, ergbeta1[idx], type=’l’, ylab=’Values of beta1’, xlab=’Iterations’,main=’(d) Ergodic Mean Plot of beta1’, ylim=ylims1)lines(idx2, ergbeta12[idx2-burnin], col=2, lty=2)apply(beta[(burnin+1):m,],2,mean)apply(beta[(burnin+1):m,],2,sd)

注意到在獨立抽樣方法產(chǎn)生的鏈中, β0和β1有很強的負相關性:

cor(beta[(burnin+1):m,1],beta[(burnin+1):m,2])=-0.954.

這是鏈收斂很慢的原因. 在高度相關的空間上使用獨立的提議分布導致鏈的混合效率低下。

M-H算法: 多元正態(tài)提議分布

在獨立抽樣中, 我們使用的提議分布是相互獨立的, 從而導致鏈的混合效率低下. 因此, 自然而然的我們 可以考慮非獨立的提議分布. 提議分布相關陣應該和后驗分布的相關陣類似. 為此, 可以考慮利用 Fisher信息陣H(β), 提議分布取為

由似然函數(shù), 容易計算得到Fisher信息陣為

算法如下：

R語言演示代碼

在計算exp(η)時, 由于在R中, exp(710)=Inf, 因此我們對η>700, 取近似log(1 + exp(η)) ≈ η.

對鏈的收斂診斷可以看出, 鏈混合的效率很高。

逐分量的M-H算法

在M-H算法中, 按分量進行逐個更新, 其優(yōu)勢在于應用方便, 不需要考慮調(diào)節(jié)參數(shù). 算法如下:

R語言演示代碼

y<-wais[,2]x<-wais[,1]m<-10000beta0<-c(0,0) #initial valuemu.beta<-c(0,0) # priors.beta<-c(100,100) # priorprop.s<-c(1.75,0.2) # sd of proposal normalbeta <- matrix(nrow=m, ncol=2)acc.prob <-c(0,0)current.beta<-beta0for (t in 1:m){    for (j in 1:2){        prop.beta<- current.beta        prop.beta[j]<- rnorm( 1, current.beta[j], prop.s[j] )        cur.eta <-current.beta[1]+current.beta[2]*x        prop.eta<-prop.beta[1]+prop.beta[2]*x        if(sum(prop.eta>700)>0) {print(t); stop;}        if(sum(cur.eta >700)>0) {print(t); stop;}        loga <-(sum(y*prop.eta-log(1+exp(prop.eta)))        -sum(y*cur.eta-log(1+exp(cur.eta)))        +sum(dnorm(prop.beta,mu.beta,s.beta,log=TRUE))        -sum(dnorm(current.beta,mu.beta,s.beta,log=TRUE)))        u<-runif(1)        u<-log(u)        if(u< loga){            current.beta<-prop.beta        acc.prob[j] <- acc.prob[j]+1        }    }    beta[t,]<-current.beta}print(acc.prob/m)

鏈的收斂診斷顯示了相比于獨立的M-H算法, 混合的效率要高的多。

III、 吉布斯采樣(Gibbs Sampling)

吉布斯采樣算法是Stuart Geman 和Donald Geman 兩兄弟在1984年提出來的，以物理學家Josiah Willard Gibbs的名字命名。這個算法在現(xiàn)代貝葉斯分析中占據(jù)重要位置。

在大數(shù)據(jù)時代，數(shù)據(jù)特征非常的多，M-H采樣由于接受概率α(i,j)的存在，在高維時需要的計算時間非常的可觀，同時α(i,j)一般小于1，有時候辛苦計算出來卻被拒絕了。能不能做到不拒絕轉(zhuǎn)移呢？吉布斯采樣算法就這樣出場了。

吉布斯采樣算法是單分量M-H算法的特殊情況，但是更容易實現(xiàn)，因而被廣泛使用。

吉布斯采樣特別適合于采樣貝葉斯網(wǎng)絡的后驗分布，因為貝葉斯網(wǎng)絡通常被指定為條件分布的集合。

吉布斯采樣算法是單分量M-H算法的特殊情況，主要有兩點特殊：

1> 吉布斯采樣算法的提議分布是當前變量xj，j = 1,2,…，k的滿條件概率分布

2> 吉布斯采樣算法的接受概率α = 1

一般的采樣方法，生成一維隨機變量并不困難，對于高維分布

生成各分量非獨立的隨機向量就非常困難。吉布斯采樣法很擅長解決這類問題，其廣泛用于多元變量聯(lián)合分布的抽樣和估計。

吉布斯采樣法把一次一個n維樣本變?yōu)椤皀次一維”樣本，但前提是我們事先知道完全條件概率分布

數(shù)學上稱為半共軛，為方便常簡寫為

不同于M-H算法，該算法不引入新的提議分布q( . )，僅通過條件分布得到以給定分布p(x)為平穩(wěn)分布的馬氏鏈的轉(zhuǎn)移概率。

把條件分布當作提議分布q( . )，則有

接受概率

以接受概率為1的方式實現(xiàn)鏈的滿概率跳轉(zhuǎn)。

吉布斯采樣構(gòu)造的馬氏鏈收斂到平穩(wěn)分布的證明

圖5-8

假設有一個二維概率分布p(x,y), 考察x坐標相同的兩個點A(x1,y1), B(x1, y2)。如圖5-8所示。

由條件概率公式推得

上面兩式右邊相等，所以

故

從上面等式，可以看出，在x=x1這條平行于y軸的直線上，如果使用條件分布p(y|x1)作為任意兩點之間的轉(zhuǎn)移概率，那么任意兩點之間的轉(zhuǎn)移滿足細致平穩(wěn)條件。同樣地，如果在y=y1這條直線上任意取兩個點A(x1,y1), C(x2, y1), 也有如下等式：

于是，我們可以構(gòu)造平面上任意兩點之間的轉(zhuǎn)移概率矩陣Q

根據(jù)上式，我們很容易驗證對平面上任意兩點X、Y (如圖5-8的兩點A、B), 滿足細致平穩(wěn)條件

因此，這個二維空間上的馬氏鏈將收斂到平穩(wěn)分布p(x,y)。

馬氏鏈的狀態(tài)轉(zhuǎn)移只是輪換地沿著x軸和y軸做轉(zhuǎn)移，于是得到樣本(x0,y0), (x0,y1), (x1,y1), (x1, y2), (x2, y2), ...馬氏鏈收斂后，最終得到的樣本就是p(x, y)的樣本。坐標軸輪換采樣是不強制要求的，也可以在x軸和y軸之間隨機地選一個坐標軸，然后按條件概率做轉(zhuǎn)移，馬氏鏈也是一樣收斂的。

我們很容易把二維推廣到高維的情形，對于二維細致平穩(wěn)條件

如果x1變?yōu)槎嗑SX1，可以看出推導過程不變，所以多維細致平穩(wěn)條件同樣是成立的。

此時轉(zhuǎn)移矩陣Q由條件分布p(y|X1)定義。上式只是說明了一個坐標軸的情況，和二維情形類似，很容易驗證對所有的坐標軸都有類似的結(jié)論。

所以，n維空間對于概率分布p(x1,x2,...,xn)可以定義如下轉(zhuǎn)移矩陣

1> 如果當前狀態(tài)為(x1,x2,...,xn)，馬氏鏈狀態(tài)轉(zhuǎn)移的過程中，只能沿著坐標軸做轉(zhuǎn)移，轉(zhuǎn)移概率定義為

2> 其它無法沿著坐標軸進行的跳轉(zhuǎn)，轉(zhuǎn)移概率都設置為0。

于是，我們可以把Gibbs Sampling算法從二維的概率分布p(x,y)推廣到n維的概率分布p(x1,x2,...,xn)。

二維Gibbs采樣算法步驟：

例子 已知：貝葉斯后驗分布有如下形式

求：利用Gibbs采樣法獲取滿足該后驗分布的樣本序列

解：

由上面兩式，求得條件概率分布

這里A(y), B(x)是關于y, x的函數(shù)，當指定y, x時，它們是常數(shù)。

求得的條件概率分布隨機變量x, y都符合正態(tài)分布，即

按已知正態(tài)分布f(x|y,data)和f(y|x,data)在兩個坐標軸上不停輪換采樣就可以得到后驗分布的樣本集。

圖5-9 Bibbs采樣點軌跡圖

從圖5-9的實驗結(jié)果可以看出采樣點大都集中在高密區(qū)。

小結(jié)：

由于MCMC的馬氏鏈的固有性質(zhì) - 采樣點依賴前一個狀態(tài)，這決定了相鄰的采樣點之間并不獨立，存在自相關。自相關的存在將降低采樣精度，為了減少樣本的依賴性，我們在實驗中可以采取如下措施：

• 刪減(Thinning)技術：對于生成的馬氏鏈，每隔一段距離保留一個，其余舍棄；

• 塊采樣(Blocked Gibbs Sampling)技術；

• 折疊采樣（Collapsed Gibbs Sampling)技術；

• 有序的過松弛(Ordered Overrelaxation)。

高維Gibbs采樣算法步驟：

當采樣空間的維數(shù)比較高時，接受-拒絕采樣和重要性采樣都不適用，原因是當維數(shù)升高時很難找到適合的提議分布，采樣效率比較差。沒有比較就沒有傷害，高維Gibbs采樣算法的優(yōu)勢就體現(xiàn)出來了，其優(yōu)點為：

• 適用于廣泛且困難的問題；

• 當聯(lián)合分布不明確或很難直接采樣，但是每個變量的條件分布是已知的并且很容易（或至少更容易）采樣時，Gibbs采樣更有效；

• 問題維數(shù)的增加通常不會降低其收斂速度或使其更復雜。

六、MCMC收斂診斷方法

MCMC方法依賴于模擬的收斂性，下面介紹三種常用的收斂性診斷方法。

1、同時產(chǎn)生多條馬爾科夫鏈

這種方法的思路是選取多個不同的初值，同時產(chǎn)生多條馬爾科夫鏈，經(jīng)過一段時間后，若這幾條鏈穩(wěn)定下來，則說明算法收斂了。在實際操作中，可在同一個二維圖中畫出這些不同的馬爾科夫鏈產(chǎn)生的后驗樣本值對迭代次數(shù)的散點圖，如果經(jīng)過若干次迭代后，這些散點圖基本穩(wěn)定，重合在一起，則可判斷其算法收斂。

2、比率診斷法

這種方法的思路是選取多個不同的初值，同時產(chǎn)生T條馬爾科夫鏈，記第j條鏈的標準差估計為Sj，鏈內(nèi)方差的均值為W，鏈間方差為B，則：

R的值趨近1，則表明MCMC模擬收斂，且比較穩(wěn)定。

3、GewekeTest法

GewekeTest由一系列的Z檢驗組成，其基本思路是：先對所有樣本(假設有M個)做一次Z檢驗，然后去掉最前面的c個樣本，用剩余的M-c個樣本重復上述檢驗，再繼續(xù)去掉最前面的c個樣本，用剩余的M-2c個樣本重復上述檢驗，依次類推，重復K次這樣的檢驗，直到M-Kc。

七、附錄

附錄1 證明n-步轉(zhuǎn)移矩陣等于其之前所有轉(zhuǎn)移矩陣的連續(xù)矩陣相乘

該證明內(nèi)容我們用定理來表述

定理： 對于?n≥0，?i，j∈S, (S是狀態(tài)空間)，則有

證明：

該定理可通過數(shù)學歸納法和Chapman-Kolmogorov等式來證明。

大家都知道，對于矩陣乘法，下面式子顯然成立。

其中n≥0, m≥0.

Chapman-Kolmogorov等式

其中n≥0, m≥0， i，j∈S，S是狀態(tài)空間。

下面只需證明Chapman-Kolmogorov等式， 然后就可以用數(shù)學歸納法證明我們的定理了。

證明方法一：

證明方法二：

根據(jù)該等式用數(shù)學歸納法證明定理略。

附錄2 WinBUGS - 用MCMC進行貝葉斯推斷的專用軟件包

WinBUGS(Bayesian inference Using Gibbs Sampling)是由英國的Imperial College和MRC聯(lián)合開發(fā)的MCMC進行貝葉斯推斷的專用軟件包。在網(wǎng)站
www.mrc-bsu.cam.ac.uk/software/bugs/可以免費下載。使用WinBUGS可以很方便地對許多常用的模型和分布進行Gibbs抽樣，使用者不需要知道參數(shù)的先驗密度或似然的精確表達式，只要設置好變量的先驗分布并對所研究的模型進行一般性的描述，就能很容易實現(xiàn)對模型的貝葉斯分析，而不需要復雜的編程。在WinBUGS中可以使用有向圖模型方式對模型進行直觀的描述，也可以直接編寫模型程序，并給出參數(shù)的 Gibbs抽樣動態(tài)圖、用smoothing方法得到的后驗分布的核密度估計圖、抽樣值的自相關圖及均數(shù)和置信區(qū)間的變換圖等，使抽樣結(jié)果更直觀、可靠。Gibbs抽樣收斂后，可以得到參數(shù)的后驗分布的平均值、標準差、95%置信區(qū)間和中位數(shù)等信息。

附錄3 MCNP - 基于蒙特卡羅方法的核物理方面的通用軟件包

MCNP(Monte Carlo N Particle Transport Code)是由美國洛斯阿拉莫斯國家實驗室(LosAlamos National Laboratory)開發(fā)的基于蒙特卡羅方法的用于計算三維復雜幾何結(jié)構(gòu)中的中子、光子、電子或者耦合中子/光子/電子輸運問題的通用軟件包，也具有計算核臨界系統(tǒng)（包括次臨界和超臨界系統(tǒng)）本征值問題的能力。

目前，MCNP以其靈活、通用的特點以及強大的功能被廣泛應用于輻射防護與射線測定、輻射屏蔽設計優(yōu)化、反應堆設計、（次）臨界裝置實驗、醫(yī)學以及檢測器設計與分析等學科領域，并得到一致認可。

附錄4 R語言，S-PLUS，MATLBA - 模擬MCMC的軟件工具

附錄5 馬爾可夫發(fā)展歷史

1906年Andrey Andreyevich Markov 引入了馬爾可夫鏈的概念

Markov, A. A. (1906). Rasprostranenie zakona bol’shih chisel na velichiny, zavisyaschie drug ot druga. Izvestiya Fiziko-matematicheskogo obschestva pri Kazanskom universitete, 15(135-156), 18

1953年梅特羅波利斯將蒙特卡洛方法引入馬爾可夫鏈

Metropolis, Nicholas; Rosenbluth, Arianna W.; Rosenbluth, Marshall N.; Teller, Augusta H.; Teller, Edward (1953). 'Equation of State Calculations by Fast Computing Machines'.

1957年馬爾可夫決策過程被提出

Bellman, R. (1957). A Markovian decision process. Journal of Mathematics and Mechanics, 679-684.

1966年Leonard E. Baum等學者提出了隱馬爾可夫模型（Hidden Markov Model，HMM）

Baum, L. E.; Petrie, T. (1966).Statistical Inference for Probabilistic Functions of Finite State Markov Chains. The Annals of Mathematical Statistics. 37 (6): 1554–1563.

1975年Baker將HMM用于語音識別

Baker, J.(1975). The DRAGON system—An overview.IEEE Transactions on Acoustics, Speech, and Signal Processing. 23: 24–29.

1980年《Markov random fields and their applications》出版，詳細描述了馬爾可夫隨機場的理論和應用

Kindermann, R., & Snell, J. L. (1980). Markov random fields and their applications (Vol. 1). American Mathematical Society.

1984年Stuart和Donald Geman兄弟描述了Gibbs抽樣算法

Geman, S.; Geman, D.(1984). Stochastic Relaxation, Gibbs Distributions, and the Bayesian Restoration of Images.IEEE Transactions on Pattern Analysis and Machine Intelligence. 6 (6): 721–741.

1988年Judea Pearl提出馬爾可夫毯的概念

Pearl, J. (2014). *Probabilistic reasoning in intelligent systems: networks of plausible inference*. Morgan Kaufmann.

1989年James D. Hamilton 1989年機應用了制轉(zhuǎn)換模型

Hamilton, J. (1989). A new approach to the economic analysis of nonstationary time series and the business cycle. Econometrica. 57 (2): 357–84.

20世紀80年代HMM 開始用于分析生物序列（DNA）

Bishop, M. and Thompson, E. (1986). Maximum Likelihood Alignment of DNA Sequences.Journal of Molecular Biology. 190 (2): 159–165.

1991年Lovejoy 研究了部分可觀測馬爾可夫決策過程（POMDP）

Lovejoy, W. S. (1991).A survey of algorithmic methods for partially observed Markov decision processes.Annals of Operations Research. 28(1):47–65.

1995年D. P. Bertsekas 和 J. N. Tsitsiklis討論了一類用于不確定條件下的控制和順序決策的動態(tài)規(guī)劃方法

Bertsekas D. P. and Tsitsiklis, J. N. (1995). Neuro-dynamic programming: an overview.Proceedings of 1995 34th IEEE Conference on Decision and Control. 1: 560-564.

1998年S.Brin和L.Page提出PageRank算法

Page, L. (1998). The pagerank citation ranking : bringing order to the web. Stanford Digital Libraries Working Paper, 9(1), 1-14.

2017年Yoshua Bengio 等研究者提出了 GibbsNet

Lamb, A. et al. (2017).GibbsNet: Iterative Adversarial Inference for Deep Graphical Models.arXiv:1712.04120.

2017年Jiaming Song, Shengjia Zhao和Stefano Ermon研究了生成對抗的訓練方法來對馬爾可夫鏈（Markov chain）的轉(zhuǎn)移算子（transition operator）進行學習

Song, J.; Zhao, S.; Ermon, S. (2017).GENERATIVE ADVERSARIAL LEARNING OF MARKOV CHAINS. ICLR.

附錄6 MCMC發(fā)展歷史

1953年6月 Metropolis等人發(fā)表了MCMC方法的開篇之作<<通過快速計算器計算狀態(tài)方程>>(Equations of state calculations by fast computing machines), 論文的主要關注點是計算積分公式(類似于貝葉斯的后驗分布)

1970年Hastings完成了MCMC方法史上里程碑的論文-Monte Carlo Sampling Methods Using Markov Chains and their Applications <<馬爾可夫鏈蒙特卡洛抽樣方法及其應用>>.解決了核物理研究中碰到的粒子分布高緯計算問題，攻克了常規(guī)蒙特卡洛方法遇到的維度瓶頸問題，將Metropolis算法一般化，并推廣為一種統(tǒng)計模擬工具，形成了Metropolis-Hastings算法。

M-H算法相對于Metropolis方法而言，看起來更像專業(yè)的統(tǒng)計模擬工具。最重要的是，M-H算法不要求“提議分布函數(shù)”必須是對稱的，從而應用起來更加靈活方便。Hastings在文章中列舉了三個例子：第一個目標分布是泊松分布，采用的提議分布是隨機游走；第二個目標分布是高斯分布，提議分布是均勻隨機游走，但此均勻分布的中心不是馬氏鏈的當前值θ， 而是-θ;第三個目標分布是多元分布，Hastings引進了Gibbs抽樣的策略，每次只更新其中一個變量，這樣的轉(zhuǎn)移矩陣同樣滿足平穩(wěn)分布。

1973年Hastings的唯一的博士生Peskun發(fā)表了題為Optimum Monte-Carlo Sampling Using Markov Chains <<最優(yōu)馬爾可夫鏈蒙特卡洛抽樣方法>>的文章，比較了Metropolis和Barker的接受概率的形式，也證明了在離散情形下Metropolis是最優(yōu)的選擇，這里的最優(yōu)性可以通過經(jīng)驗平均值的漸進方差來表示。

1984年Geman S和Geman D發(fā)表了在MCMC方法史上具有重要突破性的文章-Stochastic Relaxation,Gibbs Distributions and the Bayesian Restoration of Images. 該文基于隨機松弛(Stochastic relaxation)算法，采用Gibbs分布對圖像的貝葉斯恢復進行了研究，提出了Gibbs采樣的概念并將其引入到統(tǒng)計應用領域，Robert和Casella將此文稱為“革命的種子”，吹響了MCMC方法革命的號角。

1986年Adrian Smith 做了關于分層模型的系列學術演講。

1987年Tanner和Wong在論文中采用基于多個條件分布進行模擬的方法，這種思路等價于從聯(lián)合分布進行模擬，基本具備了Gibbs采樣的雛形。

1989年Adrian Smith第一次詳細闡釋了Gibbs抽樣的本質(zhì)。

1990年Gelfand和Smith發(fā)表論文Sampling-based approaches to calculation marginal densities.<<基于抽樣的邊際分布計算方法>> ，將Gibbs抽樣的本質(zhì)闡述得更為深刻和完整，成為主流統(tǒng)計學界大規(guī)模使用MCMC方法的真正起點。

1995年Green提出得逆跳的馬爾可夫鏈蒙特卡洛模擬(Reversible Jump Markov Chain Monte Carlo, RJMCMC)，可被視為MCMC方法第二次革命的開端，所構(gòu)建的馬氏鏈不僅可以在一個模型的參數(shù)空間內(nèi)進行轉(zhuǎn)移，還可以在不同模型(維度可以不同)之間跳躍，從而為貝葉斯模型選擇提供了強大工具。

二十世紀九十年代及以后 Richardson和Green將RJMCMC應用到混合階估計中；Brooks，Giudici和Roberts提出了提高RJMCMC轉(zhuǎn)移效率的方法；Marin和Robert將RJMCMC應用到變量選擇中。由于MCMC方法得到的并非獨立樣本，具有自相關性，因此不能直接用于參數(shù)估計，尤其是不能直接用于直接估計參數(shù)的標準誤差。為了克服自相關性，Hobert等人提出在原始馬氏鏈轉(zhuǎn)移鏈基礎上間隔抽取構(gòu)成新的樣本序列，從而克服自相關性，可以得到參數(shù)估計以及估計誤差。

八、參考文獻

[1] N. Metropolis et al., 'Equations of state calculations by fast computing machines',《the Journal of Chemical Physics》,21,1087(1953).

[2] www.random.org

[3] en.wikipedia.org/wiki/Gibbs_sampling