[Bayes] MCMC (Markov Chain Monte Carlo)

不錯的文章：LDA-math-MCMC 和 Gibbs Sampling

可做爲精進MCMC抽樣方法的學習材料。

 

 

簡單機率分佈的模擬

Box-Muller變換原理詳解

本質上來講，計算機只能生產符合均勻分佈的採樣。若是要生成其餘分佈的採樣，就須要藉助一些技巧性的方法，例如咱們在前面的文章提到過的逆變換採樣、拒絕採樣以及自適應的拒絕採樣等等。

涉及到 "逆變換" [Bayes] runif: Inversion Sampling

例如：U1, U2是均勻分佈，可獲得兩個高斯分佈的變量X, Y。

 

 

 

複雜機率分佈的模擬

使用的必要性

當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算法就是最經常使用的一種，這兩個方法在現代貝葉斯分析中被普遍使用。要了解這兩個算法，咱們首先要對馬氏鏈的平穩分佈的性質有基本的認識。

 

馬氏鏈及其平穩分佈

平穩性：這個收斂行爲主要是由機率轉移矩陣P決定的。

 

天然的，這個收斂現象並不是是咱們這個馬氏鏈獨有的，而是絕大多數馬氏鏈的共同行爲，關於馬氏鏈的收斂咱們有以下漂亮的定理：

馬氏鏈定理： 若是一個非週期馬氏鏈具備轉移機率矩陣P,且它的任何兩個狀態是連通的，那麼 limn→∞Pnij 存在且與i無關，記 limn→∞Pnij=π(j), 咱們有

1. 1. limn→∞Pn=⎡⎣⎢⎢⎢⎢⎢π(1)π(1)⋯π(1)⋯π(2)π(2)⋯π(2)⋯⋯⋯⋯⋯⋯π(j)π(j)⋯π(j)⋯⋯⋯⋯⋯⋯⎤⎦⎥⎥⎥⎥⎥
   2.  π(j)=∑i=0∞π(i)Pij
   3. π 是方程 πP=π 的惟一非負解

其中,  π=[π(1),π(2),⋯,π(j),⋯],∑i=0∞πi=1

π稱爲馬氏鏈的平穩分佈。 

這個馬氏鏈的收斂定理很是重要，全部的 MCMC(Markov Chain Monte Carlo) 方法都是以這個定理做爲理論基礎的。 

 

歷史由來 

馬氏鏈的平穩分佈 --> Metropolis算法

對於給定的機率分佈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算法也被遴選爲二十世紀的十個最重要的算法之一。

 

改進變種：Metropolis-Hastings 算法

咱們接下來介紹的MCMC 算法是 Metropolis 算法的一個改進變種，即經常使用的 Metropolis-Hastings 算法。

 

Gibbs Sampling

對於，因爲接受率 α的存在(一般 α<1), 以上 Metropolis-Hastings 算法的效率不夠高。可否找到一個轉移矩陣Q使得接受率 α=1 呢？