应用随机过程张波商豪_Markov链的应用一：MCMC算法

本文是张迪同学对马尔链的应用的介绍

应用一：Markov链在MCMC算法中的应用

1. MCMC概念

MCMC即马尔科夫链蒙特卡洛方法(Markov Chain Monte Carlo)。该方法将Markov过程引入到Monte Carlo模拟中，实现抽样分布随模拟的进行而改变的动态模拟，弥补了传统的蒙特卡罗积分只能静态模拟的缺陷。

1.1 Monte Carlo方法

基本原理是：通过设计某种实验，得出某个特定事件发生的频率，用这个频率来近似表示这一事件发生的概率，从而得到问题的数值解。例如在求积分

时，如果我们很难求解出f(x)的原函数，那么这个积分比较难求解。当然我们可以通过蒙特卡罗方法来模拟求解近似值。最简单的方法是在[a,b]之间随机的采样一个点。比如x₀，然后用f(x₀)代表在[a,b]区间上所有的f(x)的值。那么上面的定积分的近似求解为: (b−a)f(x₀)。
这种方法虽然可以求解出近似的解，但是它隐含了一个假定，即x在[a,b]之间是均匀分布的，因此在绝大部分情况，用上面的方法，模拟求出的结果很可能和真实值相差甚远。
怎么解决这个问题呢？如果我们可以得到x在[a,b]的概率分布函数p(x)，那么我们的定积分求和可以这样进行：上式最右边的这个形式就是蒙特卡罗方法的一般形式。当然这里是连续函数形式的蒙特卡罗方法，但是在离散时一样成立。可以看出，最上面我们假设x在[a,b]之间是均匀分布的时候，p(x_i)=1/(b−a)，带入我们有概率分布的蒙特卡罗积分的上式，可以得到：那么问题就变成了如何求出x的分布p(x)对应的若干个样本上来，而Markov链就是帮助找到这些复杂概率分布的对应的采样样本集的得力助手，特别是在维数非常高的情况。

1.2 Markov链的部分性质

1.2.1Markov链定义：

随机过程{,n=0,1,2……}称为Markov链，若它只取有限或可列个值，并且对任意n≥0，及任意状态i，j，，……，有

上式刻画了Markov链的特性，即“马氏性(无后效性)”，即将来的状态只依赖于现在的状态，与过去的状态独立。

1.2.2Markov链部分性质：

除马氏性外，马尔可夫链可能具有不可约性、常返性、周期性和遍历性。一个不可约和正常返的马尔可夫链是严格平稳的马尔可夫链，拥有唯一的平稳分布。遍历马尔可夫链的极限分布收敛于其平稳分布。

1.3 MCMC方法步骤

MCMC是一种对高维随机向量抽样的方法，即模拟一个马氏链，使马氏链的平稳分布为目标分布，由此产生大量的近似服从目标分布的样本，但样本不是相互独立的。MCMC的目标分布密度函数或概率函数可以只计算到差一个常数倍的值。
MCMC方法的理论依据是几个极限定理：
首先,由《应用随机过程》书中定理5.3.3可知,遍历的Markov链(即不可约、正常返、非周期的Markov链)的极限分布是平稳分布且是唯一的平稳分布。设正常返的不可约马氏链的平稳分布为, 设h()是状态空间上的有界函数，则

这类似于独立同分布随机变量平均值的强大数律。当Y服从平稳分布时，上式中的极限就等于Eh(Y)。所以，若能设计转移概率矩阵满足正常返、不可约性质的马氏链且其平稳分布为,则从某个初始分布出发按照转移概率模拟一列马氏链，由上式可以估计关于Y~的随机变量的函数的期望Eh(Y)。
因此，MCMC方法首先应建立一个Markov链,使得其极限分布是平稳分布,那么从目标分布中产生随机样本,就是从达到平稳状态的Markov链中产生样本路径。
这种方法的关键在于如何从第t时刻转移到第t+1时刻。好的转移算法应该使得马氏链比较快地收敛到平稳分布,并且不会长时间地停留在取值空间的局部区域内。
最后总结一下MCMC的采样过程：

1)输入我们任意选定的马尔科夫链状态转移矩阵Q，平稳分布π(x)，设定状态转移次数阈值，需要的样本个数；

2)从任意简单概率分布采样得到初始状态值

3)for t=0 to n1+n2−1:

a) 从条件概率分布Q(x|)中采样得到样本

b) 从均匀分布采样u∼uniform[0,1]

c) 如果u,)=π()Q(,)，则接受转移→，即=

d) 否则不接受转移，即=样本集(,,...,),即为我们需要的平稳分布对应的样本集。

这个算法虽然理论上可行，但比较难在实际中应用，因为第三步的c步骤的接受率α(,)可能非常的小，比如0.1，导致我们大部分的采样值都被拒绝转移，采样效率很低。有可能我们采样了上百万次马尔可夫链还没有收敛，也就是上面这个n1要非常非常的大。这就需要MCMC改进版采样: M-H采样和Gibbs采样来解决实际问题了。

2. M-H抽样和GiBBS抽样

M-H采样是Metropolis-Hastings采样的简称，这个算法首先由Metropolis提出，被Hastings改进，因此被称之为Metropolis-Hastings采样或M-H采样,M-H采样解决了我们上一节MCMC采样接受率过低的问题。
我们回到MCMC采样的 [细致平稳条件](即如果非周期马尔科夫链的状态转移矩阵P和概率分布对于所有的i,j满足：,则称概率分布是状态转移矩阵P的平稳分布。具体证明略。)

我们采样效率低的原因是α(i,j)太小了，比如为0.1，而α(j,i)为0.2。即：
这时我们可以看到，如果两边同时扩大五倍，接受率提高到了0.5，但是细致平稳条件却仍然是满足的，即：这样我们的接受率可以做如下改进，即：
通过这个微小的改造，我们就得到了可以在实际应用中使用的M-H采样算法过程如下：

1)输入我们任意选定的马尔科夫链状态转移矩阵Q，平稳分布π(x)，设定状态转移次数阈值n1，需要的样本个数n2

2)从任意简单概率分布采样得到初始状态值

3)for t=0 to n1+n2−1:

a) 从条件概率分布Q(x|)中采样得到样本

b) 从均匀分布采样u∼uniform[0,1]

c) 如果u,)=min(ππ,1),则接受转移→，即=

d) 否则不接受转移，即=
样本集(xn1,xn1+1,...,xn1+n2−1)即为我们需要的平稳分布对应的样本集。
很多时候，我们选择的马尔科夫链状态转移矩阵Q如果是对称的，即满足Q(i,j)=Q(j,i),这时我们的接受率可以进一步简化为：

αππ
下面我们将用一个简单的例子(普通的一维正态分布其实没必要用这种方法，仅供示意)去展示M-H采样。
我们的目标平稳分布是一个均值3，标准差2的正态分布，而选择的Markov链状态转移矩阵Q(i,j)的条件转移概率是以i为均值,方差1的正态分布在位置j的值。代码如下：

import randomimport mathfrom scipy.stats import normimport matplotlib.pyplot as plt%matplotlib inline

def norm_dist_prob(theta):    y = norm.pdf(theta, loc=3, scale=2)    return y

T = 5000pi = [0 for i in range(T)]sigma = 1t = 0while t -1:    t = t + 1    pi_star = norm.rvs(loc=pi[t - 1], scale=sigma, size=1, random_state=None)    alpha = min(1, (norm_dist_prob(pi_star[0]) / norm_dist_prob(pi[t - 1])))

    u = random.uniform(0, 1)    if u         pi[t] = pi_star[0]    else:        pi[t] = pi[t - 1]

plt.scatter(pi, norm.pdf(pi, loc=3, scale=2))num_bins = 50plt.hist(pi, num_bins, density=True, facecolor='red', alpha=0.7)plt.show()

输出结果见下图，可以看到采样值的分布与真实的分布之间的关系，采样集还是比较拟合对应分布的。

M-H采样完整解决了使用MCMC方法需要的任意概率分布样本集的问题，因此在实际生产环境得到了广泛的应用。但是在大数据时代，M-H采样面临着两大难题：

1) 我们的数据特征非常的多，M-H采样由于接受率计算式ππ的存在，在高维时需要的计算时间非常的可观，算法效率很低。同时α(i,j)一般小于1，有时候辛苦计算出来却被拒绝了。能不能做到不拒绝转移呢？

2) 由于特征维度大，很多时候我们甚至很难求出目标的各特征维度联合分布，但是可以方便求出各个特征之间的条件概率分布。这时候我们能不能只有各维度之间条件概率分布的情况下方便的采样呢？
Gibbs采样解决了上面两个问题，因此在当下大数据时代，MCMC采样基本主要是实施Gibbs采样，限于篇幅，如果有同学感兴趣的话，可以参考《应用随机过程》教材第十一章来进一步了解。

7.参考文献

a).统计计算.Dongfeng LI
b).MCMC(一)蒙特卡罗方法，马尔科夫链，MCMC采样和M-H采样.刘建平Pinard博客
c).《应用随机过程》。张波，商豪.