取樣方法（二）MCMC相關演算法介紹及程式碼實現

0.引子

書接前文，在取樣方法（一）中我們講到了拒絕取樣、重要性取樣一系列的蒙特卡洛取樣方法，但這些方法在高維空間時都會遇到一些問題，因為很難找到非常合適的可採樣Q分佈，同時保證取樣效率以及精準度。
本文將會介紹取樣方法中最重要的一族演算法，MCMC（Markov Chain Monte Carlo），在之前我們的蒙特卡洛模擬都是按照如下公式進行的：

E[f(x)]≈1m∑i=1mf(xi).xi∼p.iid
我們的x都是獨立取樣出來的，而在MCMC中，它變成了
E[f(x)]≈1m∑i=1mf(xi).(x0,x1,...,xm)∼MC(p)
其中的MC(p)就是我們本文的主角之一，馬爾可夫過程

（Markov Process）生成的馬爾可夫鏈（Markov Chain）。

1.Markov Chain（馬爾可夫鏈）

在序列的演算法（一·a）馬爾可夫模型中我們就說到了馬爾可夫模型的馬爾可夫鏈，簡單來說就是滿足馬爾可夫假設

P(sn|s0,s1,...,sn−1)=P(sn|sn−1)
的狀態序列，由馬爾可夫過程（Markov Process）生成。
一個馬爾可夫過程由兩部分組成，一是狀態集合Ω,二是轉移概率矩陣T。
其流程是：選擇一個初始的狀態分佈π0，然後進行狀態的轉移：
πn=πn−1T
得到的π0,π1,π2...πn狀態分佈序列。

注意：我們在這裡講的理論和推導都是基於離散變數的，但其實是可以直接推廣到連續變數

。

馬爾可夫鏈在很多場景都有應用，比如大名鼎鼎的pagerank演算法，都用到了類似的轉移過程；
馬爾可夫鏈在某種特定情況下，有一個奇妙的性質：

在某種條件下，當你從一個分佈π0開始進行概率轉移，無論你一開始π0的選擇是什麼，最終會收斂到一個固定的分佈π，叫做穩態（steady-state）。
穩態滿足條件：

π=πT
這裡可以參考《LDA數學八卦0.4.2》的例子，非常生動地描述了社會階層轉化的一個例子，也對MCMC作了非常好的講解

書歸正傳，回到我們取樣的場景，我們知道，取樣的難點就在於概率密度函式過於複雜而無法進行有效取樣，如果我們可以設計一個馬爾可夫過程，使得它最終收斂的分佈是我們想要取樣的概率分佈，不就可以解決我們的問題了麼。

前面提到了在某種特定情況下，這就是所有MCMC演算法的理論基礎Ergodic Theorem：
如果一個離散馬爾可夫鏈(x0,x1...xm)是一個與時間無關的Irreducible的離散，並且有一個穩態分佈π,則：

E[f(x)]≈1m∑i=1mf(xi).x∼π
它需要滿足的條件有這樣幾個，我們直接列在這裡，不作證明：

1.Time homogeneous: 狀態轉移與時間無關，這個很好理解。
2.Stationary Distribution: 最終是會收斂到穩定狀態的。
3.Irreducible: 任意兩個狀態之間都是可以互相到達的。
4.Aperiodic：馬爾可夫序列是非週期的，我們所見的絕大多數序列都是非週期的。

雖然這裡要求是離散的馬爾可夫鏈，但實際上對於連續的場景也是適用的，只是轉移概率矩陣變成了轉移概率函式。

2.MCMC

在上面馬爾可夫鏈中我們的所說的狀態都是某個可選的變數值，比如社會等級上、中、下，而在取樣的場景中,特別是多元概率分佈，並不是量從某個維度轉移到另一個維度，比如一個二元分佈，二維平面上的每一個點都是一個狀態，所有狀態的概率和為1! 這裡比較容易產生混淆，一定小心。

在這裡我們再介紹一個概念：
Detail balance：一個馬爾可夫過程是細緻平穩的，即對任意a,b兩個狀態：

π(a)Tab=π(b)Tba細緻平穩條件也可以推匯出一個非週期的馬爾可夫鏈是平穩的，因為每次轉移狀態i從狀態j獲得的量與j從i獲得的量是一樣的，那毫無疑問最後πT=π.

所以我們的目標就是需要構造這樣一個馬爾可夫鏈，使得它最後能夠收斂到我們期望的分佈π，而我們的狀態集合其實是固定的，所以最終目標就是構造一個合適的T，就大功告成了。

一般來說我們有:

π(x)=π̃ (x)Z
其中Z是歸一化引數Z=Σx′π̃ (x′)，因為我們通常能夠很方便地計算分子π̃ (x)，但是分母的計算因為要列舉所有的狀態所以過於複雜而無法計算。我們希望最終取樣出來的樣本符合π分佈。

2.1.Metropolis

2.1.1原理描述

Metropolis演算法算是MCMC的開山鼻祖了，它這裡假設我們已經有了一個狀態轉移概率函式T來表示轉移概率，T(a,b)表示從a轉移到b的概率(這裡T的選擇我們稍後再說),顯然通常情況下一個T是不滿足細緻平穩條件的：

π(a)T(a,b)≠π(b)T(b,a)
所以我們需要進行一些改造，加入一項Q使得等式成立：
π(a)T(a,b)Q(a,b)=π(b)T(b,a)Q(b,a)
基於對稱的原則，我們直接讓
Q(a,b)=π(

取樣方法（二）MCMC相關演算法介紹及程式碼實現

0.引子

1.Markov Chain（馬爾可夫鏈）

2.MCMC

2.1.Metropolis

2.1.1原理描述

取樣方法（二）MCMC相關演算法介紹及程式碼實現

【百度】大型網站的HTTPS實踐（二）——HTTPS加密演算法介紹

（二）k-means演算法原理以及python實現

篩素數方法（二）—— 費馬小定理及MR素數判斷

最近最久未使用（LRU）頁面置換演算法原理及模擬實現

spring原始碼剖析（二）Spring預設標籤解析及註冊實現

劃分方法聚類（二）K-MEANS演算法的改進

python中關於操作時間的方法（二）：使用datetime模塊

關於Android滑動沖突的解決方法（二）

JavaScript 字符串方法（二）

Jmeter腳本錄制方法（二）——手工編寫腳本（jmeter與fiddler結合使用）

測試人員必學的軟件快速測試方法（二）

面向對象進階——內置方法（二）

排序方法（二）補

電平轉換方法（二）

真正掌握vuex的使用方法（二）現在有的小夥伴是不是在想，以後如果我要在這裏寫自己的計算屬性怎麽辦？怎麽辦？咱們可以通過對象合並的方法去實現。通過Object.assign()合並對象：

css垂直居中幾種方法（二）

java後臺面試題整理及解答（二）JVM相關篇

spring-boot-2.0.3不一樣系列之源碼篇 - run方法（二）之prepareEnvir

JDBC資料庫連線池連線資料庫及資料庫操作DAO層設計通用更新及查詢方法（二）

取樣方法（二）MCMC相關演算法介紹及程式碼實現

0.引子

1.Markov Chain（馬爾可夫鏈）

2.MCMC

2.1.Metropolis

2.1.1原理描述

相關推薦