對於一般的分布的采樣，在很多的編程語言中都有實現，如最基本的滿足均勻分布的隨機數，但是對於復雜的分布，要想對其采樣，卻沒有實現好的函數，在這里，可以使用馬爾可夫鏈蒙特卡羅(Markov Chain Monte Carlo, MCMC)方法，其中Metropolis-Hasting采樣和Gibbs采樣是MCMC中使用較為廣泛的兩種形式。
MCMC的基礎理論為馬爾可夫過程，在MCMC算法中，為了在一個指定的分布上采樣，根據馬爾可夫過程，首先從任一狀態出發，模擬馬爾可夫過程，不斷進行狀態轉移，最終收斂到平穩分布。
一、馬爾可夫鏈
1、馬爾可夫鏈
設$X_t$表示隨機變量$X$在離散時間$t$時刻的取值。若該變量隨時間變化的轉移概率僅僅依賴於它的當前取值，即
$$P\left ( X_{t+1}=s_j\mid X_0=s_{0},X_1=s_{1}, \cdots ,X_t=s_i \right )=P\left ( X_{t+1}=s_j\mid X_t=s_i \right )$$
也就是說狀態轉移的概率只依賴於前一個狀態。稱這個變量為馬爾可夫變量，其中，$s_0,s_1,\cdots ,s_i,s_j\in \Omega $為隨機變量$X$可能的狀態。這個性質稱為馬爾可夫性質，具有馬爾可夫性質的隨機過程稱為馬爾可夫過程。
馬爾可夫鏈指的是在一段時間內隨機變量$X$的取值序列$\left ( X_0,X_1,\cdots ,X_m \right )$，它們滿足如上的馬爾可夫性質。
2、轉移概率
馬爾可夫鏈是通過對應的轉移概率定義的，轉移概率指的是隨機變量從一個時刻到下一個時刻，從狀態$s_i$轉移到另一個狀態$s_j$的概率，即：
$$P\left ( i\rightarrow j \right ):=P_{i,j}=P\left ( X_{t+1}=s_j\mid X_t=s_i \right )$$
記$\pi _k^{\left ( t \right )}$表示隨機變量$X$在時刻$t$的取值為$s_k$的概率，則隨機變量$X$在時刻$t+1$的取值為$s_i$的概率為：
$$\begin{align}
\pi i^{\left ( t+1 \right )} &=P\left ( X{t+1}=s_i \right ) \
&= \sum_{k}P\left ( X_{t+1}=s_i\mid X_{t}=s_k \right )\cdot P\left ( X_{t}=s_k \right )\
&= \sum_{k}P_{k,i}\cdot \pi _k^{\left ( t \right )}
\end{align}$$
假設狀態的數目為$n$，則有：
$$\left ( \pi 1^{\left ( t+1 \right )},\cdots ,\pi n^{\left ( t+1 \right )} \right )=\left ( \pi 1^{\left ( t \right )},\cdots ,\pi n^{\left ( t \right )} \right )\begin{bmatrix}
P{1,1} & P{1,2} & \cdots & P{1,n}\
P{2,1} & P_{2,2} & \cdots & P_{2,n}\
\vdots & \vdots & & \vdots \
P_{n,1} & P_{n,2} & \cdots & P_{n,n}
\end{bmatrix}$$
3、馬爾可夫鏈的平穩分布
對於馬爾可夫鏈，需要注意以下的兩點：
- 1、周期性：即經過有限次的狀態轉移，又回到了自身；
- 2、不可約：即兩個狀態之間相互轉移；
如果一個馬爾可夫過程既沒有周期性，又不可約，則稱為各態遍歷的。
對於一個各態遍歷的馬爾可夫過程，無論初始值$\pi ^{\left ( 0 \right )}$取何值，隨着轉移次數的增多，隨機變量的取值分布最終都會收斂到唯一的平穩分布$\pi ^{\ast }$，即：
$$\underset{t\rightarrow \infty }{lim}\pi ^{\left ( 0 \right )}\mathbf{P}^t=\pi ^{\ast }$$
且這個平穩分布$\pi ^{\ast }$滿足：
$$\pi ^{\ast }\mathbf{P}=\pi ^{\ast }$$
其中，$\mathbf{P}=\left ( p_{i,j} \right )_{n\times n}$為轉移概率矩陣。
二、馬爾可夫鏈蒙特卡羅方法
1、基本思想
對於一個給定的概率分布$P\left (X \right )$，若是要得到其樣本，通過上述的馬爾可夫鏈的概念，我們可以構造一個轉移矩陣為$\mathbf{P}$的馬爾可夫鏈，使得該馬爾可夫鏈的平穩分布為$P\left (X \right )$，這樣，無論其初始狀態為何值，假設記為$x_0$，那么隨着馬爾科夫過程的轉移，得到了一系列的狀態值，如：$x_0,x_1,x_2,\cdots ,x_n,x_{n+1},\cdots ,$，如果這個馬爾可夫過程在第$n$步時已經收斂，那么分布$P\left (X \right )$的樣本即為$x_n,x_{n+1},\cdots $。
2、細致平穩條件
對於一個各態遍歷的馬爾可夫過程，若其轉移矩陣為$\mathbf{P}$，分布為$\pi \left ( x \right )$，若滿足：
$$\pi \left ( i \right )P_{i,j}=\pi \left ( j \right )P_{j,i}$$
則$\pi \left ( x \right )$是馬爾可夫鏈的平穩分布，上式稱為細致平穩條件。
3、Metropolis采樣算法
Metropolis采樣算法是最基本的基於MCMC的采樣算法。
3.1、Metropolis采樣算法的基本原理
假設需要從目標概率密度函數$p\left ( \theta \right )$中進行采樣，同時，$\theta $滿足$-\infty <\theta <\infty $。Metropolis采樣算法根據馬爾可夫鏈去生成一個序列：
$$\theta ^{\left ( 1 \right )}\rightarrow \theta ^{\left ( 2 \right )}\rightarrow \cdots \theta ^{\left (t \right )}\rightarrow $$
其中，$ \theta ^{\left (t \right )}$表示的是馬爾可夫鏈在第$t$代時的狀態。
在Metropolis采樣算法的過程中，首先初始化狀態值$\theta ^{\left (1 \right )}$，然后利用一個已知的分布$q\left ( \theta \mid \theta ^{\left ( t-1 \right )} \right )$生成一個新的候選狀態$\theta ^{\left (\ast \right )}$，隨后根據一定的概率選擇接受這個新值，或者拒絕這個新值，在Metropolis采樣算法中，概率為：
$$\alpha =min: \left ( 1,; \frac{p\left ( \theta ^{\left ( \ast \right )} \right )}{p\left ( \theta ^{\left ( t-1 \right )} \right )} \right )$$
這樣的過程一直持續到采樣過程的收斂，當收斂以后，樣本$\theta ^{\left (t \right )}$即為目標分布$p\left ( \theta \right )$中的樣本。
3.2、Metropolis采樣算法的流程
基於以上的分析，可以總結出如下的Metropolis采樣算法的流程：
- 初始化時間$t=1$
- 設置$u$的值，並初始化初始狀態$\theta ^{\left (t \right )}=u$
- 重復一下的過程：
- 令$t=t+1$
- 從已知分布$q\left ( \theta \mid \theta ^{\left ( t-1 \right )} \right )$中生成一個候選狀態$\theta ^{\left (\ast \right )}$
- 計算接受的概率：$\alpha =min: \left ( 1,; \frac{p\left ( \theta ^{\left ( \ast \right )} \right )}{p\left ( \theta ^{\left ( t-1 \right )} \right )} \right )$
- 從均勻分布$Uniform\left ( 0, 1 \right )$生成一個隨機值$a$
- 如果$a\leqslant \alpha $，接受新生成的值：$\theta ^{\left (t \right )}=\theta ^{\left (\ast \right )}$；否則：$\theta ^{\left (t \right )}=\theta ^{\left (t-1 \right )}$
- 直到$t=T$
3.3、Metropolis算法的解釋
要證明Metropolis采樣算法的正確性，最重要的是要證明構造的馬爾可夫過程滿足如上的細致平穩條件，即：
$$\pi \left ( i \right )P_{i,j}=\pi \left ( j \right )P_{j,i}$$
對於上面所述的過程，分布為$p\left ( \theta \right )$，從狀態$i$轉移到狀態$j$的轉移概率為：
$$P_{i,j} =\alpha {i,j}\cdot Q{i,j}$$
其中，$Q_{i,j}$為上述已知的分布。
對於選擇該已知的分布，在Metropolis采樣算法中，要求該已知的分布必須是對稱的，即$Q_{i,j}=Q_{j,i}$，即
$$q\left ( \theta =\theta ^{\left ( t \right )}\mid \theta ^{\left ( t-1 \right )} \right )=q\left ( \theta =\theta ^{\left ( t-1 \right )}\mid \theta ^{\left ( t \right )} \right )$$
常用的符合對稱的分布主要有：正態分布，柯西分布以及均勻分布等。
接下來，需要證明在Metropolis采樣算法中構造的馬爾可夫鏈滿足細致平穩條件。
$$\begin{align}
p\left ( \theta ^{\left ( i \right )} \right )P_{i,j} &=p\left ( \theta ^{\left ( i \right )} \right )\cdot \alpha {i,j}\cdot Q{i,j} \
&= p\left ( \theta ^{\left ( i \right )} \right )\cdot min; \left { 1,\frac{p\left ( \theta ^{\left ( j \right )} \right )}{p\left ( \theta ^{\left ( i \right )} \right )} \right }\cdot Q_{i,j}\
&=min; \left { p\left ( \theta ^{\left ( i \right )} \right )Q_{i,j},p\left ( \theta ^{\left ( j \right )} \right )Q_{i,j} \right }\
&=p\left ( \theta ^{\left ( j \right )} \right )\cdot min; \left { \frac{p\left ( \theta ^{\left ( i \right )} \right )}{p\left ( \theta ^{\left ( j \right )} \right )}, 1 \right }\cdot Q_{j,i}\
&=p\left ( \theta ^{\left ( j \right )} \right )\cdot \alpha {j,i}\cdot Q{j,i}\
&=p\left ( \theta ^{\left ( j \right )} \right )P_{j,i}
\end{align}$$
因此，通過以上的方法構造出來的馬爾可夫鏈是滿足細致平穩條件的。
3.4、實驗
假設需要從柯西分布中采樣數據，我們利用Metropolis采樣算法來生成樣本，其中，柯西分布的概率密度函數為：
$$f\left ( \theta \right )=\frac{1}{\pi \left ( 1+\theta ^2 \right )}$$
那么，根據上述的Metropolis采樣算法的流程，接受概率$\alpha $的值為：
$$\alpha =min; \left ( 1,\frac{1+\left [ \theta ^{\left ( t \right )} \right ]^2}{1+\left [ \theta ^{\left ( \ast \right )} \right ]^2} \right )$$
代碼如下：
參考文獻
- 1、馬爾可夫鏈蒙特卡羅算法
- 2、受限玻爾茲曼機（RBM）學習筆記（一）預備知識
- 3、LDA數學八卦