1、隨機(jī)模擬

隨機(jī)模擬方法有一個(gè)很酷的別名是蒙特卡羅方法贾节。這個(gè)方法的發(fā)展始于20世紀(jì)40年代瘩将。
統(tǒng)計(jì)模擬中有一個(gè)很重要的問(wèn)題就是給定一個(gè)概率分布p(x),我們?nèi)绾卧谟?jì)算機(jī)中生成它的樣本鞋拟，一般而言均勻分布的樣本是相對(duì)容易生成的，通過(guò)線性同余發(fā)生器可以生成偽隨機(jī)數(shù)楚殿，我們用確定性算法生成[0,1]之間的偽隨機(jī)數(shù)序列后，這些序列的各種統(tǒng)計(jì)指標(biāo)和均勻分布Uniform(0,1)的理論計(jì)算結(jié)果非常接近竿痰，這樣的偽隨機(jī)序列就有比較好的統(tǒng)計(jì)性質(zhì)脆粥，可以被當(dāng)成真實(shí)的隨機(jī)數(shù)使用。
而我們常見(jiàn)的概率分布影涉，無(wú)論是連續(xù)的還是離散的分布变隔，都可以基于Uniform(0, 1) 的樣本生成，比如正態(tài)分布可以通過(guò)著名的 Box-Muller變換得到常潮。其他幾個(gè)著名的連續(xù)分布弟胀，包括指數(shù)分布，Gamma分布，t分布等孵户，都可以通過(guò)類似的數(shù)學(xué)變換得到萧朝，不過(guò)我們并不是總這么幸運(yùn)的，當(dāng)p(x)的形式很復(fù)雜夏哭，或者p(x)是個(gè)高維分布的時(shí)候检柬，樣本的生成就可能很困難了，此時(shí)需要一些更加復(fù)雜的隨機(jī)模擬方法來(lái)生成樣本竖配，比如MCMC方法和Gibbs采樣方法何址，不過(guò)在了解這些方法之前，我們需要首先了解一下馬爾可夫鏈及其平穩(wěn)分布进胯。

2用爪、馬爾可夫鏈

馬爾可夫鏈通俗說(shuō)就是根據(jù)一個(gè)轉(zhuǎn)移概率矩陣去轉(zhuǎn)移的隨機(jī)過(guò)程（馬爾可夫過(guò)程），該隨機(jī)過(guò)程在PageRank算法中也有使用胁镐，如下圖所示：

通俗解釋的話偎血，這里的每個(gè)圓環(huán)代表一個(gè)島嶼，比如i到j(luò)的概率是pij盯漂，每個(gè)節(jié)點(diǎn)的出度概率之和=1颇玷，現(xiàn)在假設(shè)要根據(jù)這個(gè)圖去轉(zhuǎn)移，首先我們要把這個(gè)圖翻譯成如下的矩陣：

上面的矩陣就是狀態(tài)轉(zhuǎn)移矩陣就缆，我身處的位置用一個(gè)向量表示π=(i,k,j,l)假設(shè)我第一次的位置位于i島嶼帖渠，即π0=(1,0,0,0)，第一次轉(zhuǎn)移竭宰，我們用π0乘上狀態(tài)轉(zhuǎn)移矩陣P空郊，也就是π1 = π0 * P = [pii,pij,pik,pil],也就是說(shuō)，我們有pii的可能性留在原來(lái)的島嶼i羞延，有pij的可能性到達(dá)島嶼j...第二次轉(zhuǎn)移是渣淳，以第一次的位置為基礎(chǔ)的到π2 = π1 * P，依次類推下去伴箩。

有那么一種情況入愧，我的位置向量在若干次轉(zhuǎn)移后達(dá)到了一個(gè)穩(wěn)定的狀態(tài)，再轉(zhuǎn)移π向量也不變化了嗤谚，這個(gè)狀態(tài)稱之為平穩(wěn)分布狀態(tài)π*(stationary distribution)棺蛛，這個(gè)情況需要滿足一個(gè)重要的條件，就是Detailed Balance巩步。

那么什么是Detailed Balance呢旁赊？
假設(shè)我們構(gòu)造如下的轉(zhuǎn)移矩陣：
再假設(shè)我們的初始向量為π0=(1,0,0)，轉(zhuǎn)移1000次以后達(dá)到了平穩(wěn)狀態(tài)(0.625,0.3125,0.0625)椅野。
所謂的Detailed Balance就是终畅，在平穩(wěn)狀態(tài)中：

為什么滿足了Detailed Balance條件之后杖狼，我們的馬爾可夫鏈就會(huì)收斂呢？下面的式子給出了答案：

下一個(gè)狀態(tài)是j的概率妖爷，等于從各個(gè)狀態(tài)轉(zhuǎn)移到j(luò)的概率之和蝶涩，在經(jīng)過(guò)Detailed Balance條件變換之后，我們發(fā)現(xiàn)下一個(gè)狀態(tài)是j剛好等于當(dāng)前狀態(tài)是j的概率絮识，所以馬爾可夫鏈就收斂了绿聘。

3、Markov Chain Monte Carlo

對(duì)于給定的概率分布p(x),我們希望能有便捷的方式生成它對(duì)應(yīng)的樣本次舌，由于馬爾可夫鏈能夠收斂到平穩(wěn)分布熄攘，于是一個(gè)很漂亮的想法是：如果我們能構(gòu)造一個(gè)轉(zhuǎn)移矩陣偽P的馬爾可夫鏈，使得該馬爾可夫鏈的平穩(wěn)分布恰好是p(x),那么我們從任何一個(gè)初始狀態(tài)x0出發(fā)沿著馬爾可夫鏈轉(zhuǎn)移垃它，得到一個(gè)轉(zhuǎn)移序列x0鲜屏，x1烹看，x2国拇，....xn,xn+1,如果馬爾可夫鏈在第n步已經(jīng)收斂了，于是我們就得到了p(x)的樣本xn惯殊，xn+1....

好了酱吝，有了這樣的思想，我們?cè)趺床拍軜?gòu)造一個(gè)轉(zhuǎn)移矩陣土思，使得馬爾可夫鏈最終能收斂即平穩(wěn)分布恰好是我們想要的分布p(x)呢务热？我們主要使用的還是我們的細(xì)致平穩(wěn)條件（Detailed Balance），再來(lái)回顧一下：

假設(shè)我們已經(jīng)又一個(gè)轉(zhuǎn)移矩陣為Q的馬爾可夫鏈(q(i,j)表示從狀態(tài)i轉(zhuǎn)移到狀態(tài)j的概率),顯然通常情況下：

也就是細(xì)致平穩(wěn)條件不成立己儒，所以p(x)不太可能是這個(gè)馬爾可夫鏈的平穩(wěn)分布崎岂，我們可否對(duì)馬爾可夫鏈做一個(gè)改造，使得細(xì)致平穩(wěn)條件成立呢闪湾？比如我們引入一個(gè)α(i,j)冲甘，從而使得：

那么問(wèn)題又來(lái)了，取什么樣的α(i,j)可以使上等式成立呢途样？最簡(jiǎn)單的江醇，按照對(duì)稱性：

在改造Q的過(guò)程中引入的α(i,j)稱為接受率裆站，物理意義可以理解為在原來(lái)的馬爾可夫鏈上条辟，從狀態(tài)i以q(i,j)的概率跳轉(zhuǎn)到狀態(tài)j的時(shí)候黔夭，我們以α(i,j)的概率接受這個(gè)轉(zhuǎn)移，于是得到新的馬爾可夫鏈Q(jìng)'的轉(zhuǎn)移概率q(i,j)α(i,j)羽嫡。

假設(shè)我們已經(jīng)又一個(gè)轉(zhuǎn)移矩陣Q纠修，對(duì)應(yīng)的元素為q(i,j)，把上面的過(guò)程整理一下厂僧，我們就得到了如下的用于采樣概率分布p(x)的算法：

以上的MCMC算法已經(jīng)做了很漂亮的工作了扣草，不過(guò)它有一個(gè)小問(wèn)題，馬爾可夫鏈Q(jìng)在轉(zhuǎn)移的過(guò)程中接受率α(i,j)可能偏小颜屠，這樣采樣的話容易在原地踏步辰妙，拒絕大量的跳轉(zhuǎn)，這是的馬爾可夫鏈便利所有的狀態(tài)空間要花費(fèi)太長(zhǎng)的時(shí)間甫窟，收斂到平穩(wěn)分布p(x)的速度太慢密浑，有沒(méi)有辦法提升一些接受率呢？當(dāng)然有辦法粗井，把α(i,j)和α(j,i)同比例放大尔破，不打破細(xì)致平穩(wěn)條件就好了呀，但是我們又不能無(wú)限的放大浇衬，我們可以使得上面兩個(gè)數(shù)中最大的一個(gè)放大到1懒构，這樣我們就提高了采樣中的跳轉(zhuǎn)接受率，我們?nèi)耘擂。?/p>

于是經(jīng)過(guò)這么微小的改造胆剧，我們就得到了Metropolis-Hastings算法，該算法的步驟如下：

4醉冤、Gibbs采樣

對(duì)于高維的情形秩霍，由于接受率的存在，Metropolis-Hastings算法的效率不夠高蚁阳，能否找到一個(gè)轉(zhuǎn)移矩陣Q使得接受率α=1呢铃绒？我們從二維的情形入手，假設(shè)有一個(gè)概率分布p(x,y)螺捐，考察x坐標(biāo)相同的兩個(gè)點(diǎn)A(x1,y1) ,B(x1,y2)颠悬，我們發(fā)現(xiàn)：

基于以上等式，我們發(fā)現(xiàn)归粉，在x=x1這條平行于y軸的直線上椿疗，如果使用條件分布p(y|x1)作為任何兩個(gè)點(diǎn)之間的轉(zhuǎn)移概率，那么任何兩個(gè)點(diǎn)之間的轉(zhuǎn)移滿足細(xì)致平穩(wěn)條件糠悼，同樣的届榄，在y=y1這條平行于x軸的直線上，如果使用條件分布p(x|y1) 作為倔喂，那么任何兩個(gè)點(diǎn)之間的轉(zhuǎn)移也滿足細(xì)致平穩(wěn)條件铝条。于是我們可以構(gòu)造平面上任意兩點(diǎn)之間的轉(zhuǎn)移概率矩陣Q：

有了上面的轉(zhuǎn)移矩陣Q靖苇，我們很容易驗(yàn)證對(duì)平面上任意兩點(diǎn)X，Y班缰，滿足細(xì)致平穩(wěn)條件：

由二維的情形我們很容易推廣到高維的情形：

所以高維空間中的GIbbs 采樣算法如下：