MCMC方法总结

大家好，欢迎来到IT知识分享网。

摘 要 作为一种随机采样方法，马尔科夫链蒙特卡罗（Markov Chain Monte Carlo，以下简称MCMC）在机器学习，深度学习以及自然语言处理等领域都有广泛的应用，是很多复杂算法求解的基础。比如分解机（Factorization Machines）推荐算法，还有前面讲到的受限玻尔兹曼机（RBM）原理总结，都用到了MCMC来做一些复杂运算的近似求解。本文将对MCMC的原理做一个总结。
关键词 蒙特卡罗方法；马尔科夫链；MCMC采样；M-H采样；Gibbs采样

1 蒙特卡罗方法

1.1 MCMC概述
从名字我们可以看出，MCMC由MC组成，即蒙特卡罗方法（Monte Carlo Simulation，简称MC）和马尔科夫链（Markov Chain，简称MC）。要弄懂MCMC的原理我们首先得搞清楚蒙特卡罗方法和马尔可夫链的原理。 将用三大部分来完整学习MCMC，在本章将关注蒙特卡罗方法。
1.2 蒙特卡罗方法引入
蒙特卡罗原来是一个赌场的名称，用它作为名字大概是因为蒙特卡罗方法是一种随机模拟的方法，这很像赌博场里面的扔骰子的过程。最早的蒙特卡罗方法都是为了求解一些不太好求解的求和或者积分问题。比如积分：

如果我们很难求解出f(x)的原函数，那么这个积分比较难求解。当然我们可以通过蒙特卡罗方法来模拟求解近似值。如何模拟呢？假设我们函数图像如下图：

当然，用一个值代表[a,b]区间上所有的f(x)的值，这个假设太粗糙。那么我们可以采样[a,b]区间的n个值：
，用它们的均值来代表[a,b]区间上所有的f(x)的值。这样我们上面的定积分的近似求解为：

2 马尔科夫链

2.1 马尔可夫链概述
马尔可夫链定义本身比较简单，它假设某一时刻状态转移的概率只依赖于它的前一个状态。举个形象的比喻，假如每天的天气是一个状态的话，那么今天是不是晴天只依赖于昨天的天气，而和前天的天气没有任何关系。当然这么说可能有些武断，但是这样做可以大大简化模型的复杂度，因此马尔科夫链在很多时间序列模型中得到广泛的应用，比如循环神经网络RNN(Recurrent Neural Network)，隐式马尔可夫模型HMM(Hidden Markov Model)等，当然MCMC也需要它。
如果用精确的数学定义来描述，则假设我们的序列状态是
，那么我们在时刻 的状态的条件概率仅仅依赖于时刻 ，即：

既然某一时刻状态转移的概率只依赖于它的前一个状态，那么我们只要能求出系统中任意两个状态之间的转换概率，这个马尔科夫链的模型就定了。我们来看下图这个马尔科夫链模型的具体例子（源于维基百科）。

import numpy as np matrix= np.matrix([[0.9,0.075,0.025],[0.15,0.8,0.05],[0.25,0.25,0.5]], dtype=float) vector1= np.matrix([[0.3,0.4,0.3]], dtype=float) for i in range(100): vector1 = vector1*matrix print "Current round:" , i+1 print vector1 

部分输出结果如下：

Current round: 1 [[ 0.405 0.4175 0.1775]] Current round: 2 [[ 0.4715 0.40875 0.11975]] Current round: 3 [[ 0.5156 0.3923 0.0921]] Current round: 4 [[ 0.54591 0. 0.078555]] …… Current round: 58 [[0. 0. 0.0625 ]] Current round: 59 [[ 0. 0.3125 0.0625 ]] Current round: 60 [[ 0.625 0.3125 0.0625]] …… Current round: 99 [[ 0.625 0.3125 0.0625]] Current round: 100 [[ 0.625 0.3125 0.0625]] 

matrix= np.matrix([[0.9,0.075,0.025],[0.15,0.8,0.05],[0.25,0.25,0.5]], dtype=float) vector1= np.matrix([[0.7,0.1,0.2]], dtype=float) for i in range(100): vector1 = vector1*matrix print "Current round:" , i+1 print vector1 

部分输出结果如下：

Current round: 1 [[ 0.695 0.1825 0.1225]] Current round: 2 [[ 0.6835 0.22875 0.08775]] Current round: 3 [[ 0.6714 0.2562 0.0724]] Current round: 4 [[ 0.66079 0. 0.065795]] …… Current round: 55 [[ 0. 0. 0.0625 ]] Current round: 56 [[ 0. 0. 0.0625 ]] Current round: 57 [[ 0.625 0.3125 0.0625]] …… Current round: 99 [[ 0.625 0.3125 0.0625]] Current round: 100 [[ 0.625 0.3125 0.0625]] 

matrix = np.matrix([[0.9,0.075,0.025],[0.15,0.8,0.05],[0.25,0.25,0.5]], dtype=float) for i in range(10): matrix = matrix*matrix print "Current round:" , i+1 print matrix 

输出结果如下：

Current round: 1 [[ 0.8275 0.13375 0.03875] [ 0.2675 0.66375 0.06875] [ 0.3875 0.34375 0.26875]] Current round: 2 [[ 0.73555 0. 0.051675] [ 0.42555 0. 0.074475] [ 0.51675 0. 0.]] …… Current round: 5 [[ 0. 0. 0.0] [ 0. 0. 0.0] [ 0. 0. 0.0]] Current round: 6 [[ 0.625 0.3125 0.0625] [ 0.625 0.3125 0.0625] [ 0.625 0.3125 0.0625]] 

3 MCMC采样和M-H采样

3.1 马尔可夫链的细致平稳条件
在解决从平稳分布 ，找到对应的马尔科夫链状态转移矩阵P之前，我们还需要先看看马尔科夫链的细致平稳条件。定义如下：
如果非周期马尔科夫链的状态转移矩阵P和概率分布 对于所有的i，j满足：

1. 输入任意选定的马尔科夫链状态转移矩阵Q，平稳分布 ，设定状态转移次数阈值 ，需要的样本个数
2. 从任意简单概率分布采样得到初始状态值 ；
   
   

3.4 M-H采样实例
为了更容易理解，这里给出一个M-H采样的实例。
在例子里，我们的目标平稳分布是一个均值3，标准差2的正态分布，而选择的马尔可夫链状态转移矩阵 的条件转移概率是以i为均值,方差1的正态分布在位置j的值。这个例子仅仅用来让大家加深对M-H采样过程的理解。毕竟一个普通的一维正态分布用不着去用M-H采样来获得样本。
代码如下：

import random import math from scipy.stats import norm import matplotlib.pyplot as plt %matplotlib inline def norm_dist_prob(theta): y = norm.pdf(theta, loc=3, scale=2) return y T = 5000 pi = [0 for i in range(T)] sigma = 1 t = 0 while t < 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 < alpha: pi[t] = pi_star[0] else: pi[t] = pi[t - 1] plt.scatter(pi, norm.pdf(pi, loc=3, scale=2)) num_bins = 50 plt.hist(pi, num_bins, normed=1, facecolor='red', alpha=0.7) plt.show() 

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

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

1. 我们的数据特征非常的多，M-H采样由于接受率计算式 的存在，在高维时需要的计算时间非常的可观，算法效率很低。同时 一般小于1，有时候辛苦计算出来却被拒绝了。能不能做到不拒绝转移呢？
2. 由于特征维度大，很多时候我们甚至很难求出目标的各特征维度联合分布，但是可以方便求出各个特征之间的条件概率分布。这时候我们能不能只有各维度之间条件概率分布的情况下方便的采样呢？
   Gibbs采样解决了上面两个问题，因此在大数据时代，MCMC采样基本是Gibbs采样的天下，下一篇我们就来讨论Gibbs采样。
   

4 MCMC采样和M-H采样

4.2 二维Gibbs采样
利用上一节找到的状态转移矩阵，我们就得到了二维Gibbs采样，这个采样需要两个维度之间的条件概率。具体过程如下：

用下图可以很直观的看出，采样是在两个坐标轴上不停的轮换的。当然，坐标轴轮换不是必须的，我们也可以每次随机选择一个坐标轴进行采样。不过常用的Gibbs采样的实现都是基于坐标轴轮换的。

4.3 多维Gibbs采样
上面的这个算法推广到多维的时候也是成立的。比如一个n维的概率分布 ，我们可以通过在n个坐标轴上轮换采样，来得到新的样本。对于轮换到的任意一个坐标轴 上的转移，同马尔科夫链的状态转移概率为

from mpl_toolkits.mplot3d import Axes3D from scipy.stats import multivariate_normal samplesource = multivariate_normal(mean=[5,-1], cov=[[1,1],[1,4]]) def p_ygivenx(x, m1, m2, s1, s2): return (random.normalvariate(m2 + rho * s2 / s1 * (x - m1), math.sqrt((1 - rho  2) * (s22)))) def p_xgiveny(y, m1, m2, s1, s2): return (random.normalvariate(m1 + rho * s1 / s2 * (y - m2), math.sqrt((1 - rho  2) * (s12)))) N = 5000 K = 20 x_res = [] y_res = [] z_res = [] m1 = 5 m2 = -1 s1 = 1 s2 = 2 rho = 0.5 y = m2 for i in xrange(N): for j in xrange(K): x = p_xgiveny(y, m1, m2, s1, s2) y = p_ygivenx(x, m1, m2, s1, s2) z = samplesource.pdf([x,y]) x_res.append(x) y_res.append(y) z_res.append(z) num_bins = 50 plt.hist(x_res, num_bins, normed=1, facecolor='green', alpha=0.5) plt.hist(y_res, num_bins, normed=1, facecolor='red', alpha=0.5) plt.title('Histogram') plt.show() 

输出的两个特征各自的分布如下：

然后我们看看样本集生成的二维正态分布，代码如下：

fig = plt.figure() ax = Axes3D(fig, rect=[0, 0, 1, 1], elev=30, azim=20) ax.scatter(x_res, y_res, z_res,marker='o') plt.show() 

输出的两个特征各自的分布如下：