随机采样方法整理与讲解（MCMC、Gibbs Sampling等）
原文：http://www.cnblogs.com/xbinworld/p/4266146.html?utm_source=tuicool&utm_medium=referral
本文是對參考資料中多篇關于sampling的內容進行總結+搬運，方便以后自己翻閱。其實參考資料中的資料寫的比我好，大家可以看一下！好東西多分享！PRML的第11章也是sampling，有時間后面寫到PRML的筆記中去：)
背景
隨機模擬也可以叫做蒙特卡羅模擬(Monte Carlo Simulation)。這個方法的發展始于20世紀40年代，和原子彈制造的曼哈頓計劃密切相關，當時的幾個大牛，包括烏拉姆、馮.諾依曼、費米、費曼、Nicholas Metropolis， 在美國洛斯阿拉莫斯國家實驗室研究裂變物質的中子連鎖反應的時候，開始使用統計模擬的方法,并在最早的計算機上進行編程實現。[3]
隨機模擬中有一個重要的問題就是給定一個概率分布 p ( x ) ，我們如何在計算機中生成它的樣本。一般而言均勻分布 ?U n i f o r m ( 0 , 1 ) 的樣本是相對容易生成的。 通過線性同余發生器可以生成偽隨機數，我們用確定性算法生成 [ 0 , 1 ] 之間的偽隨機數序列后，這些序列的各種統計指標和均勻分布 ?U n i f o r m ( 0 , 1 ) ?的理論計算結果非常接近。這樣的偽隨機序列就有比較好的統計性質，可以被當成真實的隨機數使用。
下面總結這么幾點：
1、蒙特卡洛數值積分
2、均勻分布，Box-Muller 變換
3、Monte Carlo principle
4、接受-拒絕抽樣（Acceptance-Rejection sampling)
5、重要性抽樣(Importance sampling)
6、馬爾科夫鏈，馬爾科夫穩態
7、MCMC——Metropolis-Hasting算法
8、MCMC——Gibbs Sampling算法
本文是對參考資料中多篇關于sampling的內容進行總結+搬運，方便以后自己翻閱。其實參考資料中的資料寫的比我好，大家可以看一下！好東西多分享！PRML的第11章也是sampling，有時間后面寫到PRML的筆記中去：)
背景
隨機模擬也可以叫做蒙特卡羅模擬(Monte Carlo Simulation)。這個方法的發展始于20世紀40年代，和原子彈制造的曼哈頓計劃密切相關，當時的幾個大牛，包括烏拉姆、馮.諾依曼、費米、費曼、Nicholas Metropolis， 在美國洛斯阿拉莫斯國家實驗室研究裂變物質的中子連鎖反應的時候，開始使用統計模擬的方法,并在最早的計算機上進行編程實現。[3]
隨機模擬中有一個重要的問題就是給定一個概率分布?p?(?x?)?，我們如何在計算機中生成它的樣本。一般而言均勻分布??U?n?i?f?o?r?m?(?0?,?1?)?的樣本是相對容易生成的。 通過線性同余發生器可以生成偽隨機數，我們用確定性算法生成?[0?,?1?]?之間的偽隨機數序列后，這些序列的各種統計指標和均勻分布??U?n?i?f?o?rm?(?0?,?1?)??的理論計算結果非常接近。這樣的偽隨機序列就有比較好的統計性質，可以被當成真實的隨機數使用。
下面總結這么幾點：
1、蒙特卡洛數值積分
2、均勻分布，Box-Muller 變換
3、Monte Carlo principle
4、接受-拒絕抽樣（Acceptance-Rejection sampling)
5、重要性抽樣(Importance sampling)
6、馬爾科夫鏈，馬爾科夫穩態
7、MCMC——Metropolis-Hasting算法
8、MCMC——Gibbs Sampling算法
1、蒙特卡洛數值積分
如果我們要求f(x)的積分，如
而f(x)的形式比較復雜積分不好求，則可以通過數值解法來求近似的結果。常用的方法是蒙特卡洛積分：
這樣把q(x)看做是x在區間內的概率分布，而把前面的分數部門看做一個函數，然后在q(x)下抽取n個樣本，當n足夠大時，可以用采用均值來近似：
因此只要q（x）比較容易采到數據樣本就行了。隨機模擬方法的核心就是如何對一個概率分布得到樣本，即抽樣（sampling）。下面我們將介紹常用的抽樣方法。
2、均勻分布，Box-Muller 變換
在計算機中生成[0,1]之間的偽隨機數序列，就可以看成是一種均勻分布。而隨機數生成方法有很多，最簡單的如：
當然計算機產生的隨機數都是偽隨機數，不過一般也就夠用了。
[Box-Muller 變換]?如果隨機變量 U1,U2 獨立且U1,U2～Uniform[0,1]，
則 Z0,Z1 獨立且服從標準正態分布。
3、Monte Carlo principle
Monte Carlo 抽樣計算隨即變量的期望值是接下來內容的重點：X 表示隨即變量，服從概率分布 p(x), 那么要計算 f(x) 的期望，只需要我們不停從 p(x) 中抽樣xi，然后對這些f（xi）取平均即可近似f(x)的期望。
4、接受-拒絕抽樣（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 值。這兩個問題會導致拒絕率很高，無用計算增加。
5、重要性抽樣(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！
6、馬爾科夫鏈，馬爾科夫穩態
在講蒙特卡洛方法之前，必須要先講一下馬爾科夫鏈；馬氏鏈的數學定義：
也就是說前一個狀態只與當前狀態有關，而與其他狀態無關，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)的算法。
8、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
1、蒙特卡洛數值積分
如果我們要求f(x)的積分，如
而f(x)的形式比較復雜積分不好求，則可以通過數值解法來求近似的結果。常用的方法是蒙特卡洛積分：
mathtex
這樣把q(x)看做是x在區間內的概率分布，而把前面的分數部門看做一個函數，然后在q(x)下抽取n個樣本，當n足夠大時，可以用采用均值來近似：
因此只要q（x）比較容易采到數據樣本就行了。隨機模擬方法的核心就是如何對一個概率分布得到樣本，即抽樣（sampling）。下面我們將介紹常用的抽樣方法。
2、均勻分布，Box-Muller 變換
在計算機中生成[0,1]之間的偽隨機數序列，就可以看成是一種均勻分布。而隨機數生成方法有很多，最簡單的如：
當然計算機產生的隨機數都是偽隨機數，不過一般也就夠用了。
[Box-Muller 變換] 如果隨機變量 U1,U2 獨立且U1,U2～Uniform[0,1]，
image
則 Z0,Z1 獨立且服從標準正態分布。
3、Monte Carlo principle
Monte Carlo 抽樣計算隨即變量的期望值是接下來內容的重點：X 表示隨即變量，服從概率分布 p(x), 那么要計算 f(x) 的期望，只需要我們不停從 p(x) 中抽樣xi，然后對這些f（xi）取平均即可近似f(x)的期望。
4、接受-拒絕抽樣（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 值。這兩個問題會導致拒絕率很高，無用計算增加。
5、重要性抽樣(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！
6、馬爾科夫鏈，馬爾科夫穩態
在講蒙特卡洛方法之前，必須要先講一下馬爾科夫鏈；馬氏鏈的數學定義：
image
也就是說前一個狀態只與當前狀態有關，而與其他狀態無關，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。事實上，從父代到子代，收入階層的變化的轉移概率如下
使用矩陣的表示方式，轉移概率矩陣記為
image
image
我們發現從第7代人開始，這個分布就穩定不變了，事實上，在這個問題中，從任意初始概率分布開始都會收斂到這個上面這個穩定的結果。
image
注：要求圖是聯通的（沒有孤立點），同時不存在一個聯通的子圖是沒有對外的出邊的（就像黑洞一樣）。
這個馬氏鏈的收斂定理非常重要，所有的 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)。如何能做到這一點呢？我們主要使用如下的定理。
image
image
馬氏鏈轉移和接受概率
假設我們已經有一個轉移矩陣Q(對應元素為q(i,j)), 把以上的過程整理一下，我們就得到了如下的用于采樣概率分布p(x)的算法。
image
image
8、MCMC——Gibbs Sampling算法
image
平面上馬氏鏈轉移矩陣的構造
image
image
以上算法收斂后，得到的就是概率分布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
本文是對參考資料中多篇關于sampling的內容進行總結+搬運，方便以后自己翻閱。其實參考資料中的資料寫的比我好，大家可以看一下！好東西多分享！PRML的第11章也是sampling，有時間后面寫到PRML的筆記中去：)
背景
隨機模擬也可以叫做蒙特卡羅模擬(Monte Carlo Simulation)。這個方法的發展始于20世紀40年代，和原子彈制造的曼哈頓計劃密切相關，當時的幾個大牛，包括烏拉姆、馮.諾依曼、費米、費曼、Nicholas Metropolis， 在美國洛斯阿拉莫斯國家實驗室研究裂變物質的中子連鎖反應的時候，開始使用統計模擬的方法,并在最早的計算機上進行編程實現。[3]
隨機模擬中有一個重要的問題就是給定一個概率分布?p?(?x?)?，我們如何在計算機中生成它的樣本。一般而言均勻分布??U?n?i?f?o?r?m?(?0?,?1?)?的樣本是相對容易生成的。 通過線性同余發生器可以生成偽隨機數，我們用確定性算法生成?[0?,?1?]?之間的偽隨機數序列后，這些序列的各種統計指標和均勻分布??U?n?i?f?o?rm?(?0?,?1?)??的理論計算結果非常接近。這樣的偽隨機序列就有比較好的統計性質，可以被當成真實的隨機數使用。
下面總結這么幾點：
1、蒙特卡洛數值積分
2、均勻分布，Box-Muller 變換
3、Monte Carlo principle
4、接受-拒絕抽樣（Acceptance-Rejection sampling)
5、重要性抽樣(Importance sampling)
6、馬爾科夫鏈，馬爾科夫穩態
7、MCMC——Metropolis-Hasting算法
8、MCMC——Gibbs Sampling算法
1、蒙特卡洛數值積分
如果我們要求f(x)的積分，如
而f(x)的形式比較復雜積分不好求，則可以通過數值解法來求近似的結果。常用的方法是蒙特卡洛積分：
這樣把q(x)看做是x在區間內的概率分布，而把前面的分數部門看做一個函數，然后在q(x)下抽取n個樣本，當n足夠大時，可以用采用均值來近似：
因此只要q（x）比較容易采到數據樣本就行了。隨機模擬方法的核心就是如何對一個概率分布得到樣本，即抽樣（sampling）。下面我們將介紹常用的抽樣方法。
2、均勻分布，Box-Muller 變換
在計算機中生成[0,1]之間的偽隨機數序列，就可以看成是一種均勻分布。而隨機數生成方法有很多，最簡單的如：
當然計算機產生的隨機數都是偽隨機數，不過一般也就夠用了。
[Box-Muller 變換]?如果隨機變量 U1,U2 獨立且U1,U2～Uniform[0,1]，
則 Z0,Z1 獨立且服從標準正態分布。
3、Monte Carlo principle
Monte Carlo 抽樣計算隨即變量的期望值是接下來內容的重點：X 表示隨即變量，服從概率分布 p(x), 那么要計算 f(x) 的期望，只需要我們不停從 p(x) 中抽樣xi，然后對這些f（xi）取平均即可近似f(x)的期望。
4、接受-拒絕抽樣（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 值。這兩個問題會導致拒絕率很高，無用計算增加。
5、重要性抽樣(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！
6、馬爾科夫鏈，馬爾科夫穩態
在講蒙特卡洛方法之前，必須要先講一下馬爾科夫鏈；馬氏鏈的數學定義：
也就是說前一個狀態只與當前狀態有關，而與其他狀態無關，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)的算法。
8、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
本文是對參考資料中多篇關于sampling的內容進行總結+搬運，方便以后自己翻閱。其實參考資料中的資料寫的比我好，大家可以看一下！好東西多分享！PRML的第11章也是sampling，有時間后面寫到PRML的筆記中去：)
背景
隨機模擬也可以叫做蒙特卡羅模擬(Monte Carlo Simulation)。這個方法的發展始于20世紀40年代，和原子彈制造的曼哈頓計劃密切相關，當時的幾個大牛，包括烏拉姆、馮.諾依曼、費米、費曼、Nicholas Metropolis， 在美國洛斯阿拉莫斯國家實驗室研究裂變物質的中子連鎖反應的時候，開始使用統計模擬的方法,并在最早的計算機上進行編程實現。[3]
隨機模擬中有一個重要的問題就是給定一個概率分布?p?(?x?)?，我們如何在計算機中生成它的樣本。一般而言均勻分布??U?n?i?f?o?r?m?(?0?,?1?)?的樣本是相對容易生成的。 通過線性同余發生器可以生成偽隨機數，我們用確定性算法生成?[0?,?1?]?之間的偽隨機數序列后，這些序列的各種統計指標和均勻分布??U?n?i?f?o?rm?(?0?,?1?)??的理論計算結果非常接近。這樣的偽隨機序列就有比較好的統計性質，可以被當成真實的隨機數使用。
下面總結這么幾點：
1、蒙特卡洛數值積分
2、均勻分布，Box-Muller 變換
3、Monte Carlo principle
4、接受-拒絕抽樣（Acceptance-Rejection sampling)
5、重要性抽樣(Importance sampling)
6、馬爾科夫鏈，馬爾科夫穩態
7、MCMC——Metropolis-Hasting算法
8、MCMC——Gibbs Sampling算法
1、蒙特卡洛數值積分
如果我們要求f(x)的積分，如
而f(x)的形式比較復雜積分不好求，則可以通過數值解法來求近似的結果。常用的方法是蒙特卡洛積分：
這樣把q(x)看做是x在區間內的概率分布，而把前面的分數部門看做一個函數，然后在q(x)下抽取n個樣本，當n足夠大時，可以用采用均值來近似：
因此只要q（x）比較容易采到數據樣本就行了。隨機模擬方法的核心就是如何對一個概率分布得到樣本，即抽樣（sampling）。下面我們將介紹常用的抽樣方法。
2、均勻分布，Box-Muller 變換
在計算機中生成[0,1]之間的偽隨機數序列，就可以看成是一種均勻分布。而隨機數生成方法有很多，最簡單的如：
當然計算機產生的隨機數都是偽隨機數，不過一般也就夠用了。
[Box-Muller 變換]?如果隨機變量 U1,U2 獨立且U1,U2～Uniform[0,1]，
則 Z0,Z1 獨立且服從標準正態分布。
3、Monte Carlo principle
Monte Carlo 抽樣計算隨即變量的期望值是接下來內容的重點：X 表示隨即變量，服從概率分布 p(x), 那么要計算 f(x) 的期望，只需要我們不停從 p(x) 中抽樣xi，然后對這些f（xi）取平均即可近似f(x)的期望。
4、接受-拒絕抽樣（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 值。這兩個問題會導致拒絕率很高，無用計算增加。
5、重要性抽樣(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！
6、馬爾科夫鏈，馬爾科夫穩態
在講蒙特卡洛方法之前，必須要先講一下馬爾科夫鏈；馬氏鏈的數學定義：
也就是說前一個狀態只與當前狀態有關，而與其他狀態無關，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)的算法。
8、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
本文是對參考資料中多篇關于sampling的內容進行總結+搬運，方便以后自己翻閱。其實參考資料中的資料寫的比我好，大家可以看一下！好東西多分享！PRML的第11章也是sampling，有時間后面寫到PRML的筆記中去：)
背景
隨機模擬也可以叫做蒙特卡羅模擬(Monte Carlo Simulation)。這個方法的發展始于20世紀40年代，和原子彈制造的曼哈頓計劃密切相關，當時的幾個大牛，包括烏拉姆、馮.諾依曼、費米、費曼、Nicholas Metropolis， 在美國洛斯阿拉莫斯國家實驗室研究裂變物質的中子連鎖反應的時候，開始使用統計模擬的方法,并在最早的計算機上進行編程實現。[3]
隨機模擬中有一個重要的問題就是給定一個概率分布?p?(?x?)?，我們如何在計算機中生成它的樣本。一般而言均勻分布??U?n?i?f?o?r?m?(?0?,?1?)?的樣本是相對容易生成的。 通過線性同余發生器可以生成偽隨機數，我們用確定性算法生成?[0?,?1?]?之間的偽隨機數序列后，這些序列的各種統計指標和均勻分布??U?n?i?f?o?rm?(?0?,?1?)??的理論計算結果非常接近。這樣的偽隨機序列就有比較好的統計性質，可以被當成真實的隨機數使用。
下面總結這么幾點：
1、蒙特卡洛數值積分
2、均勻分布，Box-Muller 變換
3、Monte Carlo principle
4、接受-拒絕抽樣（Acceptance-Rejection sampling)
5、重要性抽樣(Importance sampling)
6、馬爾科夫鏈，馬爾科夫穩態
7、MCMC——Metropolis-Hasting算法
8、MCMC——Gibbs Sampling算法
1、蒙特卡洛數值積分
如果我們要求f(x)的積分，如
而f(x)的形式比較復雜積分不好求，則可以通過數值解法來求近似的結果。常用的方法是蒙特卡洛積分：
這樣把q(x)看做是x在區間內的概率分布，而把前面的分數部門看做一個函數，然后在q(x)下抽取n個樣本，當n足夠大時，可以用采用均值來近似：
因此只要q（x）比較容易采到數據樣本就行了。隨機模擬方法的核心就是如何對一個概率分布得到樣本，即抽樣（sampling）。下面我們將介紹常用的抽樣方法。
2、均勻分布，Box-Muller 變換
在計算機中生成[0,1]之間的偽隨機數序列，就可以看成是一種均勻分布。而隨機數生成方法有很多，最簡單的如：
當然計算機產生的隨機數都是偽隨機數，不過一般也就夠用了。
[Box-Muller 變換]?如果隨機變量 U1,U2 獨立且U1,U2～Uniform[0,1]，
則 Z0,Z1 獨立且服從標準正態分布。
3、Monte Carlo principle
Monte Carlo 抽樣計算隨即變量的期望值是接下來內容的重點：X 表示隨即變量，服從概率分布 p(x), 那么要計算 f(x) 的期望，只需要我們不停從 p(x) 中抽樣xi，然后對這些f（xi）取平均即可近似f(x)的期望。
4、接受-拒絕抽樣（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 值。這兩個問題會導致拒絕率很高，無用計算增加。
5、重要性抽樣(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！
6、馬爾科夫鏈，馬爾科夫穩態
在講蒙特卡洛方法之前，必須要先講一下馬爾科夫鏈；馬氏鏈的數學定義：
也就是說前一個狀態只與當前狀態有關，而與其他狀態無關，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)的算法。
8、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
本文是對參考資料中多篇關于sampling的內容進行總結+搬運，方便以后自己翻閱。其實參考資料中的資料寫的比我好，大家可以看一下！好東西多分享！PRML的第11章也是sampling，有時間后面寫到PRML的筆記中去：)背景
隨機模擬也可以叫做蒙特卡羅模擬(Monte Carlo Simulation)。這個方法的發展始于20世紀40年代，和原子彈制造的曼哈頓計劃密切相關，當時的幾個大牛，包括烏拉姆、馮.諾依曼、費米、費曼、Nicholas Metropolis， 在美國洛斯阿拉莫斯國家實驗室研究裂變物質的中子連鎖反應的時候，開始使用統計模擬的方法,并在最早的計算機上進行編程實現。[3]
隨機模擬中有一個重要的問題就是給定一個概率分布?p?(?x?)?，我們如何在計算機中生成它的樣本。一般而言均勻分布??U?n?i?f?o?r?m?(?0?,?1?)?的樣本是相對容易生成的。 通過線性同余發生器可以生成偽隨機數，我們用確定性算法生成?[0?,?1?]?之間的偽隨機數序列后，這些序列的各種統計指標和均勻分布??U?n?i?f?o?rm?(?0?,?1?)??的理論計算結果非常接近。這樣的偽隨機序列就有比較好的統計性質，可以被當成真實的隨機數使用。
下面總結這么幾點：
1、蒙特卡洛數值積分
2、均勻分布，Box-Muller 變換
3、Monte Carlo principle
4、接受-拒絕抽樣（Acceptance-Rejection sampling)
5、重要性抽樣(Importance sampling)
6、馬爾科夫鏈，馬爾科夫穩態
7、MCMC——Metropolis-Hasting算法
8、MCMC——Gibbs Sampling算法
1、蒙特卡洛數值積分
如果我們要求f(x)的積分，如
而f(x)的形式比較復雜積分不好求，則可以通過數值解法來求近似的結果。常用的方法是蒙特卡洛積分：
這樣把q(x)看做是x在區間內的概率分布，而把前面的分數部門看做一個函數，然后在q(x)下抽取n個樣本，當n足夠大時，可以用采用均值來近似：
因此只要q（x）比較容易采到數據樣本就行了。隨機模擬方法的核心就是如何對一個概率分布得到樣本，即抽樣（sampling）。下面我們將介紹常用的抽樣方法。
2、均勻分布，Box-Muller 變換
在計算機中生成[0,1]之間的偽隨機數序列，就可以看成是一種均勻分布。而隨機數生成方法有很多，最簡單的如：
當然計算機產生的隨機數都是偽隨機數，不過一般也就夠用了。
[Box-Muller 變換]?如果隨機變量 U1,U2 獨立且U1,U2～Uniform[0,1]，
則 Z0,Z1 獨立且服從標準正態分布。
3、Monte Carlo principle
Monte Carlo 抽樣計算隨即變量的期望值是接下來內容的重點：X 表示隨即變量，服從概率分布 p(x), 那么要計算 f(x) 的期望，只需要我們不停從 p(x) 中抽樣xi，然后對這些f（xi）取平均即可近似f(x)的期望。
4、接受-拒絕抽樣（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 值。這兩個問題會導致拒絕率很高，無用計算增加。
5、重要性抽樣(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！
6、馬爾科夫鏈，馬爾科夫穩態
在講蒙特卡洛方法之前，必須要先講一下馬爾科夫鏈；馬氏鏈的數學定義：
也就是說前一個狀態只與當前狀態有關，而與其他狀態無關，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)的算法。
8、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
本文是對參考資料中多篇關于sampling的內容進行總結+搬運，方便以后自己翻閱。其實參考資料中的資料寫的比我好，大家可以看一下！好東西多分享！PRML的第11章也是sampling，有時間后面寫到PRML的筆記中去：)
背景
隨機模擬也可以叫做蒙特卡羅模擬(Monte Carlo Simulation)。這個方法的發展始于20世紀40年代，和原子彈制造的曼哈頓計劃密切相關，當時的幾個大牛，包括烏拉姆、馮.諾依曼、費米、費曼、Nicholas Metropolis， 在美國洛斯阿拉莫斯國家實驗室研究裂變物質的中子連鎖反應的時候，開始使用統計模擬的方法,并在最早的計算機上進行編程實現。[3]
隨機模擬中有一個重要的問題就是給定一個概率分布?p?(?x?)?，我們如何在計算機中生成它的樣本。一般而言均勻分布??U?n?i?f?o?r?m?(?0?,?1?)?的樣本是相對容易生成的。 通過線性同余發生器可以生成偽隨機數，我們用確定性算法生成?[0?,?1?]?之間的偽隨機數序列后，這些序列的各種統計指標和均勻分布??U?n?i?f?o?rm?(?0?,?1?)??的理論計算結果非常接近。這樣的偽隨機序列就有比較好的統計性質，可以被當成真實的隨機數使用。
下面總結這么幾點：
1、蒙特卡洛數值積分
2、均勻分布，Box-Muller 變換
3、Monte Carlo principle
4、接受-拒絕抽樣（Acceptance-Rejection sampling)
5、重要性抽樣(Importance sampling)
6、馬爾科夫鏈，馬爾科夫穩態
7、MCMC——Metropolis-Hasting算法
8、MCMC——Gibbs Sampling算法
1、蒙特卡洛數值積分
如果我們要求f(x)的積分，如
而f(x)的形式比較復雜積分不好求，則可以通過數值解法來求近似的結果。常用的方法是蒙特卡洛積分：
這樣把q(x)看做是x在區間內的概率分布，而把前面的分數部門看做一個函數，然后在q(x)下抽取n個樣本，當n足夠大時，可以用采用均值來近似：
因此只要q（x）比較容易采到數據樣本就行了。隨機模擬方法的核心就是如何對一個概率分布得到樣本，即抽樣（sampling）。下面我們將介紹常用的抽樣方法。
2、均勻分布，Box-Muller 變換
在計算機中生成[0,1]之間的偽隨機數序列，就可以看成是一種均勻分布。而隨機數生成方法有很多，最簡單的如：
當然計算機產生的隨機數都是偽隨機數，不過一般也就夠用了。
[Box-Muller 變換]?如果隨機變量 U1,U2 獨立且U1,U2～Uniform[0,1]，
則 Z0,Z1 獨立且服從標準正態分布。
3、Monte Carlo principle
Monte Carlo 抽樣計算隨即變量的期望值是接下來內容的重點：X 表示隨即變量，服從概率分布 p(x), 那么要計算 f(x) 的期望，只需要我們不停從 p(x) 中抽樣xi，然后對這些f（xi）取平均即可近似f(x)的期望。
4、接受-拒絕抽樣（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 值。這兩個問題會導致拒絕率很高，無用計算增加。
5、重要性抽樣(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！
6、馬爾科夫鏈，馬爾科夫穩態
在講蒙特卡洛方法之前，必須要先講一下馬爾科夫鏈；馬氏鏈的數學定義：
也就是說前一個狀態只與當前狀態有關，而與其他狀態無關，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)的算法。
8、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
總結
以上是生活随笔為你收集整理的随机采样方法整理与讲解（MCMC、Gibbs Sampling等）的全部內容，希望文章能夠幫你解決所遇到的問題。
- 上一篇： 单高斯分布模型GSM，高斯混合模型GMM
- 下一篇： 论文笔记——Rich feature h