1. 隨機模擬
隨機模擬(或者統計模擬)方法有一個很酷的別名是蒙特卡羅方法(Monte Carlo Simulation)。這個方法的發展始於20世紀40年代，和原子彈制造的曼哈頓計划密切相關，當時的幾個大牛，包括烏拉姆、馮.諾依曼、費米、費曼、Nicholas Metropolis， 在美國洛斯阿拉莫斯國家實驗室研究裂變物質的中子連鎖反應的時候，開始使用統計模擬的方法,並在最早的計算機上進行編程實現。
隨機模擬與計算機
現代的統計模擬方法最早由數學家烏拉姆提出，被Metropolis命名為蒙特卡羅方法，蒙特卡羅是著名的賭場，賭博總是和統計密切關聯的，所以這個命名風趣而貼切，很快被大家廣泛接受。被不過據說費米之前就已經在實驗中使用了，但是沒有發表。說起蒙特卡羅方法的源頭，可以追溯到18世紀，布豐當年用於計算π的著名的投針實驗就是蒙特卡羅模擬實驗。統計采樣的方法其實數學家們很早就知道，但是在計算機出現以前，隨機數生成的成本很高，所以該方法也沒有實用價值。隨着計算機技術在二十世紀后半葉的迅猛發展，隨機模擬技術很快進入實用階段。對那些用確定算法不可行或不可能解決的問題，蒙特卡羅方法常常為人們帶來希望。
蒙特卡羅方法
統計模擬中有一個重要的問題就是給定一個概率分布p(x)，我們如何在計算機中生成它的樣本。一般而言均勻分布 Uniform(0,1)的樣本是相對容易生成的。 通過線性同余發生器可以生成偽隨機數，我們用確定性算法生成[0,1]之間的偽隨機數序列后，這些序列的各種統計指標和均勻分布 Uniform(0,1) 的理論計算結果非常接近。這樣的偽隨機序列就有比較好的統計性質，可以被當成真實的隨機數使用。
生成一個概率分布的樣本
而我們常見的概率分布，無論是連續的還是離散的分布，都可以基於Uniform(0,1)的樣本生成。例如正態分布可以通過著名的 Box-Muller 變換得到
[Box-Muller 變換] 如果隨機變量 U1,U2 獨立且U1,U2∼Uniform[0,1]，
則 Z0,Z1 獨立且服從標准正態分布。
其它幾個著名的連續分布，包括指數分布、Gamma 分布、t 分布、F 分布、Beta 分布、Dirichlet 分布等等,也都可以通過類似的數學變換得到；離散的分布通過均勻分布更加容易生成。更多的統計分布如何通過均勻分布的變換生成出來，大家可以參考統計計算的書，其中 Sheldon M. Ross 的《統計模擬》是寫得非常通俗易懂的一本。
不過我們並不是總是這么幸運的，當p(x)的形式很復雜，或者 p(x) 是個高維的分布的時候，樣本的生成就可能很困難了。 譬如有如下的情況
- p(x)=p~(x)∫p~(x)dx,而 p~(x) 我們是可以計算的，但是底下的積分式無法顯式計算。
- p(x,y) 是一個二維的分布函數，這個函數本身計算很困難，但是條件分布 p(x|y),p(y|x)的計算相對簡單;如果 p(x) 是高維的，這種情形就更加明顯。
此時就需要使用一些更加復雜的隨機模擬的方法來生成樣本。而本節中將要重點介紹的 MCMC(Markov Chain Monte Carlo) 和 Gibbs Sampling算法就是最常用的一種，這兩個方法在現代貝葉斯分析中被廣泛使用。要了解這兩個算法，我們首先要對馬氏鏈的平穩分布的性質有基本的認識。
1.2 馬氏鏈及其平穩分布
馬氏鏈的數學定義很簡單
也就是狀態轉移的概率只依賴於前一個狀態。
我們先來看馬氏鏈的一個具體的例子。社會學家經常把人按其經濟狀況分成3類：下層(lower-class)、中層(middle-class)、上層(upper-class)，我們用1,2,3 分別代表這三個階層。社會學家們發現決定一個人的收入階層的最重要的因素就是其父母的收入階層。如果一個人的收入屬於下層類別，那么他的孩子屬於下層收入的概率是 0.65, 屬於中層收入的概率是 0.28, 屬於上層收入的概率是 0.07。事實上，從父代到子代，收入階層的變化的轉移概率如下
使用矩陣的表示方式，轉移概率矩陣記為
假設當前這一代人處在下層、中層、上層的人的比例是概率分布向量 π0=[π0(1),π0(2),π0(3)]，那么他們的子女的分布比例將是 π1=π0P, 他們的孫子代的分布比例將是 π2=π1P=π0P2, ……, 第n代子孫的收入分布比例將是 πn=πn−1P=π0Pn。
假設初始概率分布為π0=[0.21,0.68,0.11]，則我們可以計算前n代人的分布狀況如下
我們發現從第7代人開始，這個分布就穩定不變了，這個是偶然的嗎？我們換一個初始概率分布π0=[0.75,0.15,0.1] 試試看，繼續計算前n代人的分布狀況如下
我們發現，到第9代人的時候, 分布又收斂了。最為奇特的是，兩次給定不同的初始概率分布，最終都收斂到概率分布 π=[0.286,0.489,0.225]，也就是說收斂的行為和初始概率分布 π0 無關。這說明這個收斂行為主要是由概率轉移矩陣P決定的。我們計算一下Pn
我們發現，當 n 足夠大的時候，這個Pn矩陣的每一行都是穩定地收斂到π=[0.286,0.489,0.225] 這個概率分布。自然的，這個收斂現象並非是我們這個馬氏鏈獨有的，而是絕大多數馬氏鏈的共同行為，關於馬氏鏈的收斂我們有如下漂亮的定理：
馬氏鏈定理： 如果一個非周期馬氏鏈具有轉移概率矩陣P,且它的任何兩個狀態是連通的，那么 limn→∞Pnij 存在且與i無關，記 limn→∞Pnij=π(j), 我們有
- limn→∞Pn=⎡⎣⎢⎢⎢⎢⎢π(1)π(1)⋯π(1)⋯π(2)π(2)⋯π(2)⋯⋯⋯⋯⋯⋯π(j)π(j)⋯π(j)⋯⋯⋯⋯⋯⋯⎤⎦⎥⎥⎥⎥⎥
- π(j)=∑i=0∞π(i)Pij
- π 是方程 πP=π 的唯一非負解
其中,
π稱為馬氏鏈的平穩分布。
這個馬氏鏈的收斂定理非常重要，所有的 MCMC(Markov Chain Monte Carlo) 方法都是以這個定理作為理論基礎的。 定理的證明相對復雜，一般的隨機過程課本中也不給證明，所以我們就不用糾結它的證明了，直接用這個定理的結論就好了。我們對這個定理的內容做一些解釋說明：
- 該定理中馬氏鏈的狀態不要求有限，可以是有無窮多個的；
- 定理中的“非周期“這個概念我們不打算解釋了，因為我們遇到的絕大多數馬氏鏈都是非周期的；
- 兩個狀態i,j是連通並非指i 可以直接一步轉移到j(Pij>0),而是指 i 可以通過有限的n步轉移到達j(Pnij>0)。馬氏鏈的任何兩個狀態是連通的含義是指存在一個n, 使得矩陣Pn 中的任何一個元素的數值都大於零。
- 我們用 Xi 表示在馬氏鏈上跳轉第i步后所處的狀態，如果 limn→∞Pnij=π(j) 存在，很容易證明以上定理的第二個結論。由於
P(Xn+1=j)=∑i=0∞P(Xn=i)P(Xn+1=j|Xn=i)=∑i=0∞P(Xn=i)Pij
上式兩邊取極限就得到 π(j)=∑i=0∞π(i)Pij
從初始概率分布 π0 出發，我們在馬氏鏈上做狀態轉移，記Xi的概率分布為πi, 則有
由馬氏鏈收斂的定理, 概率分布πi(x)將收斂到平穩分布 π(x)。假設到第n步的時候馬氏鏈收斂，則有
所以 Xn,Xn+1,Xn+2,⋯∼π(x) 都是同分布的隨機變量，當然他們並不獨立。如果我們從一個具體的初始狀態 x0 開始,沿着馬氏鏈按照概率轉移矩陣做跳轉，那么我們得到一個轉移序列 x0,x1,x2,⋯xn,xn+1⋯, 由於馬氏鏈的收斂行為， xn,xn+1,⋯都將是平穩分布 π(x) 的樣本。
1.3 Markov Chain Monte Carlo
對於給定的概率分布p(x),我們希望能有便捷的方式生成它對應的樣本。由於馬氏鏈能收斂到平穩分布， 於是一個很的漂亮想法是：如果我們能構造一個轉移矩陣為P的馬氏鏈，使得該馬氏鏈的平穩分布恰好是p(x), 那么我們從任何一個初始狀態x0出發沿着馬氏鏈轉移, 得到一個轉移序列 x0,x1,x2,⋯xn,xn+1⋯,， 如果馬氏鏈在第n步已經收斂了，於是我們就得到了 π(x) 的樣本xn,xn+1⋯。
這個絕妙的想法在1953年被 Metropolis想到了，為了研究粒子系統的平穩性質， Metropolis 考慮了物理學中常見的波爾茲曼分布的采樣問題，首次提出了基於馬氏鏈的蒙特卡羅方法，即Metropolis算法，並在最早的計算機上編程實現。Metropolis 算法是首個普適的采樣方法，並啟發了一系列 MCMC方法，所以人們把它視為隨機模擬技術騰飛的起點。 Metropolis的這篇論文被收錄在《統計學中的重大突破》中， Metropolis算法也被遴選為二十世紀的十個最重要的算法之一。
我們接下來介紹的MCMC 算法是 Metropolis 算法的一個改進變種，即常用的 Metropolis-Hastings 算法。由上一節的例子和定理我們看到了，馬氏鏈的收斂性質主要由轉移矩陣P 決定, 所以基於馬氏鏈做采樣的關鍵問題是如何構造轉移矩陣P,使得平穩分布恰好是我們要的分布p(x)。如何能做到這一點呢？我們主要使用如下的定理。
定理：[細致平穩條件] 如果非周期馬氏鏈的轉移矩陣P和分布π(x) 滿足
則 π(x) 是馬氏鏈的平穩分布，上式被稱為細致平穩條件(detailed balance condition)。
其實這個定理是顯而易見的，因為細致平穩條件的物理含義就是對於任何兩個狀態i,j, 從i 轉移出去到j 而丟失的概率質量，恰好會被從 j 轉移回i 的概率質量補充回來，所以狀態i上的概率質量π(i)是穩定的，從而π(x)是馬氏鏈的平穩分布。數學上的證明也很簡單，由細致平穩條件可得
由於π 是方程 πP=π的解，所以π是平穩分布。
假設我們已經有一個轉移矩陣為Q馬氏鏈(q(i,j)表示從狀態 i轉移到狀態j的概率，也可以寫為 q(j|i)或者q(i→j)), 顯然，通常情況下
也就是細致平穩條件不成立，所以 p(x) 不太可能是這個馬氏鏈的平穩分布。我們可否對馬氏鏈做一個改造，使得細致平穩條件成立呢？譬如，我們引入一個 α(i,j), 我們希望
取什么樣的 α(i,j) 以上等式能成立呢？最簡單的，按照對稱性，我們可以取
於是(*)式就成立了。所以有
於是我們把原來具有轉移矩陣Q的一個很普通的馬氏鏈，改造為了具有轉移矩陣Q′的馬氏鏈，而 Q′恰好滿足細致平穩條件，由此馬氏鏈Q′的平穩分布就是p(x)！
在改造 Q 的過程中引入的 α(i,j)稱為接受率，物理意義可以理解為在原來的馬氏鏈上，從狀態 i 以q(i,j) 的概率轉跳轉到狀態j 的時候，我們以α(i,j)的概率接受這個轉移，於是得到新的馬氏鏈Q′的轉移概率為q(i,j)α(i,j)。
馬氏鏈轉移和接受概率
假設我們已經有一個轉移矩陣Q(對應元素為q(i,j)), 把以上的過程整理一下，我們就得到了如下的用於采樣概率分布p(x)的算法。
上述過程中 p(x),q(x|y) 說的都是離散的情形，事實上即便這兩個分布是連續的，以上算法仍然是有效，於是就得到更一般的連續概率分布 p(x)的采樣算法，而 q(x|y) 就是任意一個連續二元概率分布對應的條件分布。
以上的 MCMC 采樣算法已經能很漂亮的工作了，不過它有一個小的問題：馬氏鏈Q在轉移的過程中的接受率 α(i,j) 可能偏小，這樣采樣過程中馬氏鏈容易原地踏步，拒絕大量的跳轉，這使得馬氏鏈遍歷所有的狀態空間要花費太長的時間，收斂到平穩分布p(x)的速度太慢。有沒有辦法提升一些接受率呢?
假設 α(i,j)=0.1,α(j,i)=0.2, 此時滿足細致平穩條件，於是
上式兩邊擴大5倍，我們改寫為
看，我們提高了接受率，而細致平穩條件並沒有打破！這啟發我們可以把細致平穩條件(**) 式中的α(i,j),α(j,i) 同比例放大，使得兩數中最大的一個放大到1，這樣我們就提高了采樣中的跳轉接受率。所以我們可以取
於是，經過對上述MCMC 采樣算法中接受率的微小改造，我們就得到了如下教科書中最常見的 Metropolis-Hastings 算法。
對於分布 p(x),我們構造轉移矩陣 Q′ 使其滿足細致平穩條件
此處 x 並不要求是一維的，對於高維空間的 p(x)，如果滿足細致平穩條件
那么以上的 Metropolis-Hastings 算法一樣有效。
1.4 Gibbs Sampling
對於高維的情形，由於接受率 α的存在(通常 α<1), 以上 Metropolis-Hastings 算法的效率不夠高。能否找到一個轉移矩陣Q使得接受率 α=1 呢？我們先看看二維的情形，假設有一個概率分布 p(x,y), 考察x坐標相同的兩個點A(x1,y1),B(x1,y2)，我們發現
所以得到
即
基於以上等式，我們發現，在 x=x1 這條平行於 y軸的直線上，如果使用條件分布 p(y|x1)做為任何兩個點之間的轉移概率，那么任何兩個點之間的轉移滿足細致平穩條件。同樣的，如果我們在 y=y1 這條直線上任意取兩個點 A(x1,y1),C(x2,y1),也有如下等式
平面上馬氏鏈轉移矩陣的構造
於是我們可以如下構造平面上任意兩點之間的轉移概率矩陣Q
有了如上的轉移矩陣 Q, 我們很容易驗證對平面上任意兩點 X,Y, 滿足細致平穩條件
於是這個二維空間上的馬氏鏈將收斂到平穩分布 p(x,y)。而這個算法就稱為 Gibbs Sampling 算法,是 Stuart Geman 和Donald Geman 這兩兄弟於1984年提出來的，之所以叫做Gibbs Sampling 是因為他們研究了Gibbs random field, 這個算法在現代貝葉斯分析中占據重要位置。
Gibbs Sampling 算法中的馬氏鏈轉移
以上采樣過程中，如圖所示，馬氏鏈的轉移只是輪換的沿着坐標軸 x軸和y軸做轉移，於是得到樣本 (x0,y0),(x0,y1),(x1,y1),(x1,y2),(x2,y2),⋯ 馬氏鏈收斂后，最終得到的樣本就是 p(x,y) 的樣本，而收斂之前的階段稱為 burn-in period。額外說明一下，我們看到教科書上的 Gibbs Sampling 算法大都是坐標軸輪換采樣的，但是這其實是不強制要求的。最一般的情形可以是，在t時刻，可以在x軸和y軸之間隨機的選一個坐標軸，然后按條件概率做轉移，馬氏鏈也是一樣收斂的。輪換兩個坐標軸只是一種方便的形式。
以上的過程我們很容易推廣到高維的情形，對於(***) 式，如果x1 變為多維情形x1，可以看出推導過程不變，所以細致平穩條件同樣是成立的
此時轉移矩陣 Q 由條件分布 p(y|x1) 定義。上式只是說明了一根坐標軸的情形，和二維情形類似，很容易驗證對所有坐標軸都有類似的結論。所以n維空間中對於概率分布 p(x1,x2,⋯,xn) 可以如下定義轉移矩陣
- 如果當前狀態為(x1,x2,⋯,xn)，馬氏鏈轉移的過程中，只能沿着坐標軸做轉移。沿着 xi 這根坐標軸做轉移的時候，轉移概率由條件概率 p(xi|x1,⋯,xi−1,xi+1,⋯,xn) 定義；
- 其它無法沿着單根坐標軸進行的跳轉，轉移概率都設置為 0。
於是我們可以把Gibbs Smapling 算法從采樣二維的 p(x,y) 推廣到采樣n 維的 p(x1,x2,⋯,xn)
以上算法收斂后，得到的就是概率分布p(x1,x2,⋯,xn)的樣本，當然這些樣本並不獨立，但是我們此處要求的是采樣得到的樣本符合給定的概率分布，並不要求獨立。同樣的，在以上算法中，坐標軸輪換采樣不是必須的，可以在坐標軸輪換中引入隨機性，這時候轉移矩陣 Q 中任何兩個點的轉移概率中就會包含坐標軸選擇的概率，而在通常的 Gibbs Sampling 算法中，坐標軸輪換是一個確定性的過程，也就是在給定時刻t，在一根固定的坐標軸上轉移的概率是1。
本文轉自：http://cos.name/2013/01/lda-math-mcmc-and-gibbs-sampling/