MCMC(一)蒙特卡羅方法
作為一種隨機采樣方法，馬爾科夫鏈蒙特卡羅（Markov Chain Monte Carlo，以下簡稱MCMC）在機器學習,深度學習以及自然語言處理等領域都有廣泛的應用，是很多復雜算法求解的基礎。比如我們前面講到的分解機(Factorization Machines)推薦算法，還有前面講到的受限玻爾茲曼機（RBM）原理總結，都用到了MCMC來做一些復雜運算的近似求解。下面我們就對MCMC的原理做一個總結。
1. MCMC概述
從名字我們可以看出，MCMC由兩個MC組成，即蒙特卡羅方法（Monte Carlo Simulation，簡稱MC）和馬爾科夫鏈（Markov Chain ，也簡稱MC）。要弄懂MCMC的原理我們首先得搞清楚蒙特卡羅方法和馬爾科夫鏈的原理。我們將用三篇來完整學習MCMC。在本篇，我們關注於蒙特卡羅方法。
2. 蒙特卡羅方法引入
蒙特卡羅原來是一個賭場的名稱，用它作為名字大概是因為蒙特卡羅方法是一種隨機模擬的方法，這很像賭博場里面的扔骰子的過程。最早的蒙特卡羅方法都是為了求解一些不太好求解的求和或者積分問題。比如積分：$$\theta = \int_a^b f(x)dx$$
如果我們很難求解出$f(x)$的原函數，那么這個積分比較難求解。當然我們可以通過蒙特卡羅方法來模擬求解近似值。如何模擬呢？假設我們函數圖像如下圖:
則一個簡單的近似求解方法是在[a,b]之間隨機的采樣一個點。比如$x_0$,然后用$f(x_0)$代表在[a,b]區間上所有的$f(x)$的值。那么上面的定積分的近似求解為:$$(b-a)f(x_0)$$
當然，用一個值代表[a,b]區間上所有的$f(x)$的值，這個假設太粗糙。那么我們可以采樣[a,b]區間的n個值：${x_0,x_1,...x_{n-1}}$,用它們的均值來代表[a,b]區間上所有的$f(x)$的值。這樣我們上面的定積分的近似求解為:$$\frac{b-a}{n}\sum\limits_{i=0}^{n-1}f(x_i)$$
雖然上面的方法可以一定程度上求解出近似的解，但是它隱含了一個假定，即$x$在[a,b]之間是均勻分布的，而絕大部分情況，$x$在[a,b]之間不是均勻分布的。如果我們用上面的方法，則模擬求出的結果很可能和真實值相差甚遠。
怎么解決這個問題呢？ 如果我們可以得到$x$在[a,b]的概率分布函數$p(x)$，那么我們的定積分求和可以這樣進行：$$\theta = \int_a^b f(x)dx = \int_a^b \frac{f(x)}{p(x)}p(x)dx \approx \frac{1}{n}\sum\limits_{i=0}^{n-1}\frac{f(x_i)}{p(x_i)}$$
上式最右邊的這個形式就是蒙特卡羅方法的一般形式。當然這里是連續函數形式的蒙特卡羅方法，但是在離散時一樣成立。
可以看出，最上面我們假設$x$在[a,b]之間是均勻分布的時候，$p(x_i) = 1/(b-a)$，帶入我們有概率分布的蒙特卡羅積分的上式，可以得到：$$\frac{1}{n}\sum\limits_{i=0}^{n-1}\frac{f(x_i)}{1/(b-a)} = \frac{b-a}{n}\sum\limits_{i=0}^{n-1}f(x_i) $$
也就是說，我們最上面的均勻分布也可以作為一般概率分布函數$p(x)$在均勻分布時候的特例。那么我們現在的問題轉到了如何求出$x$的分布$p(x)$對應的若干個樣本上來。
3. 概率分布采樣
上一節我們講到蒙特卡羅方法的關鍵是得到$x$的概率分布。如果求出了$x$的概率分布，我們可以基於概率分布去采樣基於這個概率分布的n個$x$的樣本集，帶入蒙特卡羅求和的式子即可求解。但是還有一個關鍵的問題需要解決，即如何基於概率分布去采樣基於這個概率分布的n個$x$的樣本集。
對於常見的均勻分布$uniform(0,1)$是非常容易采樣樣本的，一般通過線性同余發生器可以很方便的生成(0,1)之間的偽隨機數樣本。而其他常見的概率分布，無論是離散的分布還是連續的分布，它們的樣本都可以通過$uniform(0,1)$的樣本轉換而得。比如二維正態分布的樣本$(Z_1,Z_2)$可以通過通過獨立采樣得到的$uniform(0,1)$樣本對$(X_1,X_2)$通過如下的式子轉換而得：$$Z_1 = \sqrt{-2 ln X_1}cos(2\pi X_2)$$$$Z_2 = \sqrt{-2 ln X_1}sin(2\pi X_2)$$
其他一些常見的連續分布，比如t分布，F分布，Beta分布，Gamma分布等，都可以通過類似的方式從$uniform(0,1)$得到的采樣樣本轉化得到。在python的numpy，scikit-learn等類庫中，都有生成這些常用分布樣本的函數可以使用。
不過很多時候，我們的$x$的概率分布不是常見的分布，這意味着我們沒法方便的得到這些非常見的概率分布的樣本集。那這個問題怎么解決呢？
4. 接受-拒絕采樣
對於概率分布不是常見的分布，一個可行的辦法是采用接受-拒絕采樣來得到該分布的樣本。既然 $p(x)$ 太復雜在程序中沒法直接采樣，那么我設定一個程序可采樣的分布 $q(x)$ 比如高斯分布，然后按照一定的方法拒絕某些樣本，以達到接近 $p(x)$ 分布的目的，其中$q(x)$叫做 proposal distribution。
具體采用過程如下，設定一個方便采樣的常用概率分布函數 $q(x)$，以及一個常量 $k$，使得 $p(x)$ 總在 $kq(x)$ 的下方。如上圖。
首先，采樣得到$q(x)$的一個樣本$z_0$，采樣方法如第三節。然后，從均勻分布$（0, kq(z_0)) $中采樣得到一個值$u$。如果$u$落在了上圖中的灰色區域，則拒絕這次抽樣，否則接受這個樣本$z_0$。重復以上過程得到n個接受的樣本$z_0,z_1,...z_{n-1}$,則最后的蒙特卡羅方法求解結果為：$$\frac{1}{n}\sum\limits_{i=0}^{n-1}\frac{f(z_i)}{p(z_i)}$$
整個過程中，我們通過一系列的接受拒絕決策來達到用$q(x)$模擬$p(x)$概率分布的目的。
5. 蒙特卡羅方法小結
使用接受-拒絕采樣，我們可以解決一些概率分布不是常見的分布的時候，得到其采樣集並用蒙特卡羅方法求和的目的。但是接受-拒絕采樣也只能部分滿足我們的需求，在很多時候我們還是很難得到我們的概率分布的樣本集。比如：
1）對於一些二維分布$p(x,y)$，有時候我們只能得到條件分布$p(x|y)$和$p(y|x)$和,卻很難得到二維分布$p(x,y)$一般形式，這時我們無法用接受-拒絕采樣得到其樣本集。
2）對於一些高維的復雜非常見分布$p(x_1,x_2,...,x_n)$，我們要找到一個合適的$q(x)$和$k$非常困難。
從上面可以看出，要想將蒙特卡羅方法作為一個通用的采樣模擬求和的方法，必須解決如何方便得到各種復雜概率分布的對應的采樣樣本集的問題。而我們下一篇要講到的馬爾科夫鏈就是幫助找到這些復雜概率分布的對應的采樣樣本集的白衣騎士。下一篇我們來總結馬爾科夫鏈的原理。
（歡迎轉載，轉載請注明出處。歡迎溝通交流： liujianping-ok@163.com）