开启掘金成长之旅！这是我参与「掘金日新计划 · 2 月更文挑战」的第 32 天，点击查看活动详情

---

对于一般的分布的采样，在很多的编程语言中都有实现，如最基本的满足均匀分布的随机数，但是对于复杂的分布，要想对其采样，却没有实现好的函数，在这里，可以使用马尔可夫链蒙特卡罗(Markov Chain Monte Carlo, MCMC)方法，其中Metropolis-Hastings采样和Gibbs采样是MCMC中使用较为广泛的两种形式。

MCMC的基础理论为马尔可夫过程，在MCMC算法中，为了在一个指定的分布上采样，根据马尔可夫过程，首先从任一状态出发，模拟马尔可夫过程，不断进行状态转移，最终收敛到平稳分布。

1. 马尔可夫链

1.1. 马尔可夫链

设$X_t$表示随机变量$X$在离散时间$t$时刻的取值。若该变量随时间变化的转移概率仅仅依赖于它的当前取值，即

也就是说状态转移的概率只依赖于前一个状态。称这个变量为马尔可夫变量，其中，$s_0,s_1,\cdots ,s_i,s_j\in \Omega$为随机变量$X$可能的状态。这个性质称为马尔可夫性质，具有马尔可夫性质的随机过程称为马尔可夫过程。

马尔可夫链指的是在一段时间内随机变量$X$的取值序列$\left ( X_0,X_1,\cdots ,X_m \right )$，它们满足如上的马尔可夫性质。

1.2. 转移概率

马尔可夫链是通过对应的转移概率定义的，转移概率指的是随机变量从一个时刻到下一个时刻，从状态$s_i$转移到另一个状态$s_j$的概率，即：

记$\pi _k^{\left ( t \right )}$表示随机变量$X$在时刻$t$的取值为$s_k$的概率，则随机变量$X$在时刻$t+1$的取值为$s_i$的概率为：

假设状态的数目为$n$，则有：

1.3. 马尔可夫链的平稳分布

对于马尔可夫链，需要注意以下的两点：

• 1、周期性：即经过有限次的状态转移，又回到了自身；
• 2、不可约：即两个状态之间相互转移；

如果一个马尔可夫过程既没有周期性，又不可约，则称为各态遍历的。

对于一个各态遍历的马尔可夫过程，无论初始值$\pi ^{\left ( 0 \right )}$取何值，随着转移次数的增多，随机变量的取值分布最终都会收敛到唯一的平稳分布$\pi ^{\ast }$，即：

且这个平稳分布$\pi ^{\ast }$满足：

其中，$\mathbf{P}=\left ( p_{i,j} \right )_{n\times n}$为转移概率矩阵。

2. 马尔可夫链蒙特卡罗方法

2.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.2. 细致平稳条件

对于一个各态遍历的马尔可夫过程，若其转移矩阵为$\mathbf{P}$，分布为$\pi \left ( x \right )$，若满足：

则$\pi \left ( x \right )$是马尔可夫链的平稳分布，上式称为细致平稳条件。

2.3. Metropolis采样算法

Metropolis采样算法是最基本的基于MCMC的采样算法。

2.3.1. Metropolis采样算法的基本原理

假设需要从目标概率密度函数$p\left ( \theta \right )$中进行采样，同时，$\theta$满足$-\infty <\theta <\infty$。Metropolis采样算法根据马尔可夫链去生成一个序列：

其中，$\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采样算法中，概率为：

这样的过程一直持续到采样过程的收敛，当收敛以后，样本$\theta ^{\left (t \right )}$即为目标分布$p\left ( \theta \right )$中的样本。

2.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$

2.3.3. Metropolis算法的解释

要证明Metropolis采样算法的正确性，最重要的是要证明构造的马尔可夫过程满足如上的细致平稳条件，即：

对于上面所述的过程，分布为$p\left ( \theta \right )$，从状态$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采样算法中构造的马尔可夫链满足细致平稳条件。

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

2.3.4. 实验

假设需要从柯西分布中采样数据，我们利用Metropolis采样算法来生成样本，其中，柯西分布的概率密度函数为：

那么，根据上述的Metropolis采样算法的流程，接受概率$\alpha$的值为：

代码如下：

'''
Date:20160629
@author: zhaozhiyong
'''
import random
from scipy.stats import norm
import matplotlib.pyplot as plt

def cauchy(theta):
    y = 1.0 / (1.0 + theta ** 2)
    return y

T = 5000
sigma = 1
thetamin = -30
thetamax = 30
theta = [0.0] * (T+1)
theta[0] = random.uniform(thetamin, thetamax)

t = 0
while t < T:
    t = t + 1
    theta_star = norm.rvs(loc=theta[t - 1], scale=sigma, size=1, random_state=None)
    #print theta_star
    alpha = min(1, (cauchy(theta_star[0]) / cauchy(theta[t - 1])))
        
    u = random.uniform(0, 1)
    if u <= alpha:
        theta[t] = theta_star[0]
    else:
        theta[t] = theta[t - 1]

ax1 = plt.subplot(211)
ax2 = plt.subplot(212) 
plt.sca(ax1)
plt.ylim(thetamin, thetamax)
plt.plot(range(T+1), theta, 'g-')
plt.sca(ax2)
num_bins = 50
plt.hist(theta, num_bins, normed=1, facecolor='red', alpha=0.5)
plt.show()

实验的结果：

对于Metropolis采样算法，其要求选定的分布必须是对称的，为了弥补这样的一个缺陷，在下一篇中，介绍一下Metropolis-Hastings采样算法，其是Metropolis采样算法的推广形式。

参考文献

[1] 马尔可夫链蒙特卡罗算法

[2] 受限玻尔兹曼机（RBM）学习笔记（一）预备知识

[3] LDA数学八卦

---

开启掘金成长之旅！这是我参与「掘金日新计划 · 2 月更文挑战」的第 32 天，点击查看活动详情