简单易学的机器学习算法——马尔可夫链蒙特卡罗方法MCMC

阅读量：4967 次

发布时间：2019-06-12

本文共 5705 字，大约阅读时间需要 19 分钟。

对于一般的分布的采样，在很多的编程语言中都有实现，如最基本的满足均匀分布的随机数，但是对于复杂的分布，要想对其采样，却没有实现好的函数，在这里，可以使用马尔可夫链蒙特卡罗(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}$$

因此，通过以上的方法构造出来的马尔可夫链是满足细致平稳条件的。