如果我們要求$f(x)$的積分，可化成，
\[\int {\frac{{f(x)}}{{p(x)}}p(x)dx} \]
$p(x)$是x的概率分布，假設${g(x) = \frac{{f(x)}}{{p(x)}}}$，然后在$p(x)$的分布下，抽取x個樣本，當n足夠大時，可以采用均值來近似$f(x)$的積分，
\[\int {f(x)dx} \approx \frac{{g({x_1}) + g({x_2}) + ... + g({x_n})}}{n}\]
1. 接受-拒絕采樣
就算我們已知$p(x)$的分布，也很難得到一堆符合$p(x)$分布的樣本$\{ {x_1},{x_2},...,{x_n}\} $來帶入$g(x)$。
既然$p(x)$太復雜在程序中沒法直接采樣，那么我們設定一個程序可抽樣的分布$q(x)$比如高斯分布，然后按照一定的方法拒絕某些樣本，達到接近$p(x)$分布的目的，這就是接受拒絕采樣。
接受拒絕采樣具體操作如下，設定一個方便抽樣的函數$q(x)$，以及一個常量k，使得已知的分布$p(x)$（紅線）總在$kq(x)$（藍線）的下方，
從方便抽樣的$q(x)$分布抽樣得到$z_0$
從均勻分布$\left( {0,kq\left( {{z_0}} \right)} \right)$抽樣得到$u_0$
如果$u_0$剛好落到灰色區域，拒絕這次采樣，否則接受這次采樣$x_t=z_0$
重復以上過程，得到接近$p(x)$分布的樣本$\{ {x_1},{x_2},...,{x_n}\} $
在高維的情況下，接受-拒絕采樣會出現兩個問題，第一是合適的$q(x)$分布比較難以找到，第二是很難確定一個合理的 k 值。這兩個問題會導致拒絕率很高，無用計算增加。
2. MCMC蒙特卡洛采樣
這里我們需要引入馬爾科夫鏈來幫忙，后面的幾種采樣方法都是基於馬爾科夫鏈的，如果不想明白原理的話，這邊一段都可以跳過。
可以參照馬爾科夫細致平穩條件
我們的馬爾科夫鏈模型的狀態轉移矩陣收斂到的穩定概率分布與我們的初始狀態概率分布無關。這是一個非常好的性質，也就是說，如果我們得到了這個穩定概率分布$p(x)$對應的馬爾科夫鏈模型的狀態轉移矩陣$P$，則我們可以用任意的概率分布樣本$q_{0}(x)$開始，帶入馬爾科夫鏈模型的狀態轉移矩陣$P$，這樣經過一些序列的轉換，${q_0}\mathop \to \limits^P {q_1}\mathop \to \limits^P {q_2}...\mathop \to \limits^P {q_z}$，最終就可以得到符合對應穩定概率分布的樣本${q_z} = \{ {x_1},{x_2},...,{x_n}\} $。
馬爾科夫鏈采樣具體過程總結如下，輸入馬爾科夫狀態轉移矩陣P，需要樣本個數n，
從任意簡單概率分布得到樣本${q_0} = \{ {x_1},{x_2},...,{x_n}\}$開始
對樣本$q_t$中的每一個樣本$x_i$，利用狀態轉移矩陣P，${x_i} \times P$，求得轉移后的$q_{t+1}$
達到終止條件（轉移次數閾值或者q穩定）后的${q_z} = \{ {x_1},{x_2},...,{x_n}\} $即為我們需要的平穩分布對應的樣本集。
針對目標平穩分布$p(x)$，如何得到它所對應的馬爾科夫鏈狀態轉移矩陣P呢？
從任意一個馬爾科夫狀態轉移矩陣Q開始，一般情況，$p(x)$和Q不滿足細致平穩條件，
$\pi (i)Q(i,j) \ne \pi (j)Q(j,i)$
為了使等式成立，引入一個接受率$\alpha (i,j)$，
$\pi (i)Q(i,j)\alpha (i,j) = \pi (j)Q(j,i)\alpha (j,i)$
只要接受率$\alpha (i,j)$滿足下面兩個式子，就可以使細致平穩條件成立，
$\begin{array}{l}
\alpha (i,j) = \pi (j)Q(j,i)\\
\alpha (j,i) = \pi (i)Q(i,j)
\end{array}$
類似與接受-拒接采樣思想，我們的目標轉移矩陣P可以通過一個任意的馬爾科夫轉移矩陣Q乘以某個接受率值得到。取值在[0,1]之間，可以理解為一個概率值。接受-拒接采樣是以一個方便抽樣的分布通過一定的接受-拒絕概率得到一個不方便抽樣的分布$p(x)$樣本，MCMC采樣是以一個常見的馬爾科夫鏈狀態轉移矩陣Q通過一定的接受-拒絕概率得到目標轉移矩陣P。
目標是生成$p(x)$分布的樣本$\{ {X_1},{X_2},...,{X_n}\}$，需要樣本的個數為n，MCMC采樣具體操作如下，輸入任意選定的馬爾科夫狀態轉移矩陣Q，從任意初始值$x_0$開始，
從均勻分布采樣u∼uniform[0,1]采樣樣本u，
利用轉移矩陣Q中針對$x_{t}$元素的轉移概率，也就是條件概率分布$Q(x|x_t)$，將$x_{t}$轉移到$x_*$，
如果$u < \alpha (x_{t},x_*) = p({x_*})Q({x_*},{x_{t}})$，則接受這次轉移${x_t} \to {x_*}$，即$x_{t+1}=x_{*}$，否則不接受轉移，
重復上述步驟直至達到轉移次數閾值后，再采樣n次，得到接近$p(x)$分布的n個樣本$\{ {X_1},{X_2},...,{X_n}\}$。
這個過程基本上就是MCMC采樣的完整采樣理論了，但是這個采樣算法還是比較難在實際中應用，為什么呢？問題在於$p({x_t})Q({x_t}|{x_{t - 1}})$可能非常的小，比如0.1，導致我們大部分的采樣值都被拒絕轉移，采樣效率很低。有可能我們采樣了上百萬次馬爾可夫鏈還沒有收斂，這讓人難以接受，怎么辦呢？
3. M-H采樣
針對細致平穩條件，
$\pi (i)Q(i,j)\alpha (i,j) = \pi (j)Q(j,i)\alpha (j,i)$
接受率$\alpha (i,j)$太小怎么辦，我們等式左右兩邊同時擴大使\alpha (j,i)$達到1
$\begin{array}{l}
\alpha (i,j) = min\{ \frac{{\pi (j)Q(j,i)}}{{\pi (i)Q(i,j)}},1\} \\
\alpha (j,i) = 1
\end{array}$
M-H采樣具體操作如下，輸入任意選定的馬爾科夫狀態轉移矩陣Q，從任意初始值$x_0$開始，
從均勻分布采樣u∼uniform[0,1]采樣樣本u，
利用轉移矩陣Q中針對$x_{t}$元素的轉移概率，也就是條件概率分布$Q(x|x_t)$，將$x_{t}$轉移到$x_*$，
如果$u < \min \left\{ {\frac{{p({x_*})Q({x_*},{x_t})}}{{p({x_t})Q({x_t},{x_*})}},1} \right\}$，則接受這次轉移${x_t} \to {x_*}$，即$x_{t+1}=x_{*}$，否則不接受轉移，
重復上述步驟直至達到轉移次數閾值后，再采樣n次，得到接近$p(x)$分布的n個樣本$\{ {X_1},{X_2},...,{X_n}\}$。
import numpy as np import pandas as pd def normal(x, u, sigma): p = np.exp(-(x - u) ** 2 / (2*sigma ** 2))/(np.sqrt(2*np.pi)*sigma) return p #選擇的馬爾可夫鏈狀態轉移矩陣Q(i,j)的條件轉移概率是以i為均值,方差1的正態分布在位置j的值 #FUNCTION(x)就是我們已知函數形式，想生成該概率分布樣本的分布p(x)，可用normal(x, 0, 2)正態分布為例， # %% #從任意初始值開始，所有的樣本都存在x_t里 x_t = [0] for t in range(5000): # 從均勻分布采樣u∼uniform[0,1]采樣樣本u， u = np.random.uniform(0, 1) #利用轉移矩陣Q，將上一個時候的xt，轉移到x*（x_） x_ = np.random.normal(x_t[-1], 1) alpha = min(1, FUNCTION(x_)*normal(x_t[-1], x_, 1)) #判斷是否接受轉移 if u < alpha: x_t.append(x_) # %% # 后n個采樣樣本 即為我們需要的接近p(x)分布的n個樣本 # 查看采樣樣本的概率密度函數是否符合FUNCTION(x) pd.Series(x_t[-200:]).plot(kind='kde')
"
M-H采樣完整解決了使用蒙特卡羅方法需要的任意概率分布樣本集的問題，因此在實際生產環境得到了廣泛的應用。
但是在大數據時代，M-H采樣面臨着兩大難題：
1） 我們的數據特征非常的多，M-H采樣由於接受率計算式π(j)Q(j,i)π(i)Q(i,j)π(j)Q(j,i)π(i)Q(i,j)的存在，在高維時需要的計算時間非常的可觀，算法效率很低。同時α(i,j)α(i,j)一般小於1，有時候辛苦計算出來卻被拒絕了。能不能做到不拒絕轉移呢？
2） 由於特征維度大，很多時候我們甚至很難求出目標的各特征維度聯合分布，但是可以方便求出各個特征之間的條件概率分布。這時候我們能不能只有各維度之間條件概率分布的情況下方便的采樣呢？
Gibbs采樣解決了上面兩個問題，因此在大數據時代，MCMC采樣基本是Gibbs采樣的天下，下一篇我們就來討論Gibbs采樣。
"
4. Gibbs采樣
針對二維數據分布，
。。。。。。
由於Gibbs采樣在高維特征時的優勢，目前我們通常意義上的MCMC采樣都是用的Gibbs采樣。當然Gibbs采樣是從M-H采樣的基礎上的進化而來的，同時Gibbs采樣要求數據至少有兩個維度，一維概率分布的采樣是沒法用Gibbs采樣的,這時M-H采樣仍然成立。
有了Gibbs采樣來獲取概率分布的樣本集，有了蒙特卡羅方法來用樣本集模擬求和，他們一起就奠定了MCMC算法在大數據時代高維數據模擬求和時的作用。MCMC系列就在這里結束吧。
感恩劉建平Pinard 撒花~
http://www.cnblogs.com/pinard/p/6645766.html
論文：
可用於參數估計，
“采用馬爾科夫鏈蒙特卡洛仿真的方法對參數進行樣本抽樣，在對所抽樣的樣本進行分布擬合，即可獲得模型參數的點估計、區間估計和近似的概率分布。具體參數估計過程可以通過軟件Open BUGS 來實現。”
參考：
https://www.cnblogs.com/xbinworld/p/4266146.html
http://www.cnblogs.com/pinard/p/6625739.html
http://www.cnblogs.com/pinard/p/6645766.html