隨機採樣方法整理與講解（MCMC、Gibbs Sampling等）

本文是對參考資料中多篇關於sampling的內容進行總結+搬運，方便之後本身翻閱。其實參考資料中的資料寫的比我好，你們能夠看一下！好東西多分享！PRML的第11章也是sampling，有時間後面寫到PRML的筆記中去：)

背景

隨機模擬也能夠叫作蒙特卡羅模擬(Monte Carlo Simulation)。這個方法的發展始於20世紀40年代，和原子彈製造的曼哈頓計劃密切相關，當時的幾個大牛，包括烏拉姆、馮.諾依曼、費米、費曼、Nicholas Metropolis， 在美國洛斯阿拉莫斯國家實驗室研究裂變物質的中子連鎖反應的時候，開始使用統計模擬的方法,並在最先的計算機上進行編程實現。[3]

隨機模擬中有一個重要的問題就是給定一個機率分佈p(x)，咱們如何在計算機中生成它的樣本。通常而言均勻分佈 Uniform(0,1)的樣本是相對容易生成的。 經過線性同餘發生器能夠生成僞隨機數，咱們用肯定性算法生成[0,1]之間的僞隨機數序列後，這些序列的各類統計指標和均勻分佈 Uniform(0,1) 的理論計算結果很是接近。這樣的僞隨機序列就有比較好的統計性質，能夠被當成真實的隨機數使用。

 

下面總結這麼幾點：

一、蒙特卡洛數值積分

二、均勻分佈，Box-Muller 變換

三、Monte Carlo principle

四、接受-拒絕抽樣（Acceptance-Rejection sampling)

五、重要性抽樣(Importance sampling)

六、馬爾科夫鏈，馬爾科夫穩態

七、MCMC——Metropolis-Hasting算法

八、MCMC——Gibbs Sampling算法

 

一、蒙特卡洛數值積分

若是咱們要求f(x)的積分，如

而f(x)的形式比較複雜積分很差求，則能夠經過數值解法來求近似的結果。經常使用的方法是蒙特卡洛積分：

這樣把q(x)看作是x在區間內的機率分佈，而把前面的分數部門看作一個函數，而後在q(x)下抽取n個樣本，當n足夠大時，能夠用採用均值來近似：

所以只要q（x）比較容易採到數據樣本就好了。隨機模擬方法的核心就是如何對一個機率分佈獲得樣本，即抽樣（sampling）。下面咱們將介紹經常使用的抽樣方法。

 

二、均勻分佈，Box-Muller 變換

在計算機中生成[0,1]之間的僞隨機數序列，就能夠當作是一種均勻分佈。而隨機數生成方法有不少，最簡單的如：

固然計算機產生的隨機數都是僞隨機數，不過通常也就夠用了。

 

[Box-Muller 變換]  若是隨機變量 U1,U2 獨立且U1,U2∼Uniform[0,1]，

則 Z0,Z1 獨立且服從標準正態分佈。

 

三、Monte Carlo principle

Monte Carlo 抽樣計算隨即變量的指望值是接下來內容的重點：X 表示隨即變量，服從機率分佈 p(x), 那麼要計算 f(x) 的指望，只須要咱們不停從 p(x) 中抽樣xi，而後對這些f（xi）取平均便可近似f(x)的指望。

 

 

四、接受-拒絕抽樣（Acceptance-Rejection sampling)[2]

不少實際問題中，p(x)是很難直接採樣的的，所以，咱們須要求助其餘的手段來採樣。既然 p(x) 太複雜在程序中無法直接採樣，那麼我設定一個程序可抽樣的分佈 q(x) 好比高斯分佈，而後按照必定的方法拒絕某些樣本，達到接近 p(x) 分佈的目的，其中q(x)叫作 proposal distribution 。

具體操做以下，設定一個方便抽樣的函數 q(x)，以及一個常量 k，使得 p(x) 總在 kq(x) 的下方。（參考上圖）

• x 軸方向：從 q(x) 分佈抽樣獲得 a。(若是是高斯，就用以前說過的 tricky and faster 的算法更快）
• y 軸方向：從均勻分佈（0, kq(a)) 中抽樣獲得 u。
• 若是恰好落到灰色區域： u > p(a), 拒絕， 不然接受此次抽樣
• 重複以上過程

在高維的狀況下，Rejection Sampling 會出現兩個問題，第一是合適的 q 分佈比較難以找到，第二是很難肯定一個合理的 k 值。這兩個問題會致使拒絕率很高，無用計算增長。

 

五、重要性抽樣(Importance sampling)[2]

Importance Sampling 也是藉助了容易抽樣的分佈 q (proposal distribution)來解決這個問題，直接從公式出發：

其中，p(z) / q(z) 能夠看作 importance weight。咱們來考察一下上面的式子，p 和 f 是肯定的，咱們要肯定的是 q。要肯定一個什麼樣的分佈纔會讓採樣的效果比較好呢？直觀的感受是，樣本的方差越小指望收斂速率越快。好比一次採樣是 0, 一次採樣是 1000, 平均值是 500,這樣採樣效果不好，若是一次採樣是 499, 一次採樣是 501, 你說指望是 500,可信度還比較高。在上式中，咱們目標是 p×f/q 方差越小越好，因此 |p×f| 大的地方，proposal distribution q(z) 也應該大。舉個稍微極端的例子：

第一個圖表示 p 分佈， 第二個圖的陰影區域 f = 1，非陰影區域 f = 0, 那麼一個良好的 q 分佈應該在左邊箭頭所指的區域有很高的分佈機率，由於在其餘區域的採樣計算實際上都是無效的。這代表 Importance Sampling 有可能比用原來的 p 分佈抽樣更加有效。

可是惋惜的是，在高維空間裏找到一個這樣合適的 q 很是難。即便有 Adaptive importance sampling 和 Sampling-Importance-Resampling(SIR) 的出現，要找到一個同時知足 easy to sample 而且 good approximations 的 proposal distribution, it is often impossible！

 

六、馬爾科夫鏈，馬爾科夫穩態

在講蒙特卡洛方法以前，必需要先講一下馬爾科夫鏈；馬氏鏈的數學定義：

也就是說前一個狀態只與當前狀態有關，而與其餘狀態無關，Markov Chain 體現的是狀態空間的轉換關係，下一個狀態只決定與當前的狀態(能夠聯想網頁爬蟲原理，根據當前頁面的超連接訪問下一個網頁)。以下圖：

舉一個例子，若是當前狀態爲 u(x) = (0.5, 0.2, 0.3), 那麼下一個矩陣的狀態就是 u(x)T = (0.18, 0.64, 0.18), 依照這個轉換矩陣一直轉換下去，最後的系統就趨近於一個穩定狀態 (0.22, 0.41, 0.37) (此處只保留了兩位有效數字)。而事實證實不管你從那個點出發，通過很長的 Markov Chain 以後都會聚集到這一點。[2]

再舉一個例子，社會學家常常把人按其經濟情況分紅3類：下層(lower-class)、中層(middle-class)、上層(upper-class)，咱們用1,2,3 分別表明這三個階層。社會學家們發現決定一我的的收入階層的最重要的因素就是其父母的收入階層。若是一我的的收入屬於下層類別，那麼他的孩子屬於下層收入的機率是 0.65, 屬於中層收入的機率是 0.28, 屬於上層收入的機率是 0.07。事實上，從父代到子代，收入階層的變化的轉移機率以下

使用矩陣的表示方式，轉移機率矩陣記爲

咱們發現從第7代人開始，這個分佈就穩定不變了，事實上，在這個問題中，從任意初始機率分佈開始都會收斂到這個上面這個穩定的結果。

注：要求圖是聯通的（沒有孤立點），同時不存在一個聯通的子圖是沒有對外的出邊的（就像黑洞同樣）。

這個馬氏鏈的收斂定理很是重要，全部的 MCMC(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)。如何能作到這一點呢？咱們主要使用以下的定理。

 

 

                                    馬氏鏈轉移和接受機率

 

假設咱們已經有一個轉移矩陣Q(對應元素爲q(i,j)), 把以上的過程整理一下，咱們就獲得了以下的用於採樣機率分佈p(x)的算法。

 

 

八、MCMC——Gibbs Sampling算法

          平面上馬氏鏈轉移矩陣的構造

 

 

以上算法收斂後，獲得的就是機率分佈p(x1,x2,⋯,xn)的樣本，固然這些樣本並不獨立，可是咱們此處要求的是採樣獲得的樣本符合給定的機率分佈，並不要求獨立。一樣的，在以上算法中，座標軸輪換採樣不是必須的，能夠在座標軸輪換中引入隨機性，這時候轉移矩陣 Q 中任何兩個點的轉移機率中就會包含座標軸選擇的機率，而在一般的 Gibbs Sampling 算法中，座標軸輪換是一個肯定性的過程，也就是在給定時刻t，在一根固定的座標軸上轉移的機率是1。

 

參考資料

[1] http://blog.csdn.net/xianlingmao/article/details/7768833

[2] http://www.cnblogs.com/daniel-D/p/3388724.html

[3] http://cos.name/2013/01/lda-math-mcmc-and-gibbs-sampling/

[4] An Introduction to MCMC for Machine Learning，2003

[5] Introduction to Monte Carlo Methods