这几天在看有关multimodal learning（多模态学习）的一些内容，随之就牵扯出了boltzman机，mcmc（马尔可夫链蒙特卡洛）采样等一系列内容。总之mcmc采样是ml领域非常重要的一个理论，此外，目前能找到的书籍或网络资源很多说的都不是很清楚。所以我决定为此专门写一篇博客。
要去讲清楚一个算法或一个模型，通常首先得明确两点，一是使用它的背景，二是更具针对性的该方法能够解决的问题，或者说被设计出来的动机，即background和motivation。其次再是算法的流程与细节。本文将通过尽可能通俗易懂的语言对mcmc进行一个简单而全面的介绍。
一、background
首先，mcmc算法是用来做什么的？开门见山，mcmc就是计算机用来模拟采样的算法。那么什么是采样？采样就是给定一个分布，根据该分布从样本空间中生成样本。比如你可以通过抛硬币近似完成对p=12的0-1分布的采样。再举个最简单的通过计算机程序对低维离散样本空间进行采样的例子，比如一维变量X取值的样本空间为{1,2,3}，且取1,2,3三个值的概率分别为12,14,14。那么这时候你要如何从这个分布中进行采样？我想大多数人自己都写过这样简单的程序，首先根据各离散取值的概率大小对[0,1]区间进行等比例划分，如划分为[0,0.5],[0,5,0.75],[0.75,1]这三个区间，再通过计算机产生[0,1]之间的伪随机数，根据伪随机数的落点即可完成一次采样。好了，mcmc名字可能会让第一次见到它的人觉得很不友好，但是它其实真的只是想要单纯不做作地完成上述的采样工作。
二、motivation
前文说明了mcmc是用来采样的，也举了一个最简单的采样的例子及采样的方法。那么问题来了，如果我们要对一个连续分布（即给定一个已知的概率密度函数p(x)）进行采样，那么上述对[0,1]区间划分的方式显然就失效了（当然，最简单的方法，可能会有人考虑将概率密度函数均匀分段积分，然后继续采用之前的做法，不过这样永远只能得到近似的采样结果，永远不可能得到原始分布的采样结果，并且高维情况下积分的计算代价以及是否可积本身也是个问题）。所以，在给定概率密度函数的连续分布下要如何采样呢？mcmc是不是该登场亮相了？先别急，mcmc当然可以用来对任意的给定的概率密度函数进行采样，你甚至可以在一维离散分布时也使用mcmc来采样，所以我的意思是，割鸡焉用牛刀，暂时还不是mcmc出场的最佳时机。对于概率密度函数的采样，相信一些人可以很直观的想到可以利用累积分布函数（P(x<t)，即∫t−∞p(x)dx），即在[0,1]间随机生成一个数a，然后求使得P(x<t)=a成立的t，t即可以视作从该分部中得到的一个采样结果。但是这是基于累积分布函数可积，且P(x<t)=a可解，即累积分布函数具有反函数的条件下的。假如累积分布函数没有反函数呢？是的，这时候我们就可以使用mcmc进行采样了（当然其实还有更为直接也更为直观的方式，比如rejection sampling，不过我觉得背景介绍已经做的太多了这里不想再提了，如果感兴趣可以自行google，原理很简单）。好了，这是mcmc的motivation之一，即当累积分布函数不可积，或者不存在反函数时，mcmc可以帮助我们解决这样的抽样问题。
除此以外，考虑高维离散分布，理论上似乎你可以对所有离散分布都按文章开头所说的那种方式进行采样，但是当离散样本空间维度较高时，基于现在计算机的存储能力与计算资源，你显然不可能总是将每一个可能出现的样本与其对应的出现概率全部计算出来，并对应到[0,1]中的分段区间上，对于高维连续空间也是类似的道理。这是mcmc的另一个motivation，即旨在解决样本空间过大，计算资源过高（由其是高维样本空间组合爆炸，无论是离散还是连续）的问题。
这两点是mcmc主要的motivation，其实mcmc还被用来解决另一个问题，不过我觉得把这个问题放在本文后面的某个部分来说更加合适。
三、马尔可夫链蒙特卡洛
你现在可能对mcmc是用来解决什么问题的已经了如指掌了，但是说了半天，还是不知道mcmc到底是什么。mcmc是Markov Chain Monte Carlo的简写（我这里不打算介绍马尔可夫链或者蒙特卡洛算法，所以，如果要更舒适地阅读以下内容，最好对这二者有一定的理解。如果有时间我可能会单独写一篇马尔可夫链的博客），是的，mcmc本质上是一个马尔可夫链。那么马尔可夫链与采样有什么关系？我们很自然的会想到，马尔可夫链本身在每个时刻都有一个分布，这个分布是不是可以和我们的目标采样分布p(x)结合起来？是的，答案就是如此，我们就是要构造一个具有平稳分布(stationary distribution)p(x)（平稳分布即当时间t大于某一值时，该马尔可夫链的分布始终保持不变）的马尔可夫链，有了这样的马尔可夫链，我们要如何采样便显而易见了，我们采样的样本即马尔可夫链在达到平稳分布后时间序列上各个时刻的状态（一个状态对应一个样本）。这时候剩下的主要问题似乎就是，我们要如何构造一个具有平稳分布p(x)的马尔可夫链呢？不过我暂时先不讲这个问题，因为这里似乎还有一个问题没搞清楚，对于高维样本空间组合爆炸的问题，为什么转成用马尔可夫链就可以避免这个问题了？假如高维样本空间的样本数是n，那马尔可夫链不是要维护一个n2的状态转移矩阵？那不是更爆炸了吗？这个问题是这样的，看似每个马尔可夫链都需要一个状态转移矩阵，但是实际上我们要的是当系统处于某一状态时，在下一时刻到达其他各所能到达状态的概率。也就是说我们需要的只是一个计算当前状态xi转到下一时刻可行状态xj的概率P(xj|xi)的抽象的接口，而并不真正需要一个全局的物理的矩阵。P(xj|xi)对我们来说就像一个抽象的接口，一个输入当前状态和目标转移状态，返回概率的接口，我们只需要为当前状态和目标状态定义一个统一的转移概率计算方式就可以了。此外，我们还可以提前限定好一个状态下一时刻所能到达的可能状态数（比如高维空间状态的转移，每次转移只允许某一个维度上变化）。通过以上两种方式我们就完全避免了维度爆炸的问题。其实这里多说一句，通过mcmc来做抽样本质上就是将样本纵向在空间上的分布转化成为了样本横向在时间上的等价分布，所以mcmc具备避免空间爆炸的能力是必然的。
四、如何构建一个符合要求的马尔可夫链
现在终于可以进入mcmc采样最核心的部分，也是很多书籍或者参考资料一开头就会跟你讲的东西，即mcmc采样方法的算法本身即细节，或者直接说是构造具有与目标分布p(x)一致的平稳分布的马尔可夫链的方式。
这里首先说明一个马尔可夫链的两个重要性质，一是非周期不可约马尔可夫链一定可以收敛到唯一的平稳分布，证明可见[1]；二是马尔可夫链在某一时刻的分布π为平稳态的一个充分不必要条件是
这里πi表示处于状态xi的概率（离散情况），pij表示由状态xi转移到状态xi的概率，即P(xj|xi)。关于第二个性质，我这里简单证明一下
易见π为一个平稳分布。
所以给定目标分布p(x)要构造一个符合条件的可以用来在时间序列上进行采样的马尔可夫链，就是要构造一个马尔可夫链，使得π=p(x)且具有上述两点性质。讲到这里，其实mcmc是个什么东西已经完全讲完了，下面就是这个马尔可夫链具体如何构造的问题了，可以想象，符合条件的马尔可夫链数不胜数，接下来我主要引用两个比较经典的做法。（后面的内容，任何地方都能找到，如果不是很在意实现细节的，可以忽略）
(A) Metropolis-Hasting算法
记要采样的概率分布是π(x)，mcmc的关键是要构建出一个转移概率分布p(y|x)，使得π(x)∗p(y|x)=π(y)∗p(x|y)成立。现在假设我们有一个容易采样的分布q(y|x)（称为建议分布)，对于目前的样本x，它能够通过q(y|x)得到下一个建议样本y，这个建议样本y按照一定的概率被接受或者不被接受，称为比率α(x,y)=min{1,q(x|y)∗π(y)q(y|x)∗π(x)}。即如果知道样本xi，如何知道下一个样本xi+1是什么呢？就是通过q(y|xi)得到一个建议样本y，然后根据α(xi,y)决定xi+1=y 还是xi+1=xi。换句话说，我们这里的p(y|x)就是q(y|x)∗α(x,y)。显然，我们可以得到
这里还剩下一个问题，就是建议分布q(x)如何选取，具体参见[2].
(B) Gibbs抽样
其实Gibbs分布的思想我在第三节已经提到了，“提前限定好一个状态下一时刻所能到达的可能状态数（比如高维空间状态的转移，每次转移只允许某一个维度上变化）”这里其实说的就是一种Gibbs抽样的方法。下面我将对这种方法进行更正式地描述。样本空间为n维向量空间，目标分布仍然记为π(x)。此时，我们将状态x到状态y的转移概率定义为
这里x−i是指样本x除了第i维以外的所有维度上的取值，即x1,x2,…,xi−1,xi+1,…,xn。q(x)是样本空间为{1,2,…,n}的一个分布，一般情况下为均匀分布,V即{1,2,…,n}。此外，自身到自身的转移概率定义如下
根据上述定义，我们不难证明π(x)p(y|x)=π(y)p(x|y)（用mathjax打公式太累了，这里我就不写了）。并且，我们也容易知道该马尔可夫链是不可约的（因为每个状态都可以经过有限次或无限次转移到达任意另一状态），无周期的，所以这是一个符合要求的定义。
此外，注意到，我这里虽然限定了每次转移只有一个维度可以发生变化，但是并没有定死哪一个维度发生变化，而是通过一个概率分布来选择发生变化的维度。这样每一次采样我们还是要对每个维度的可能变化都进行采样，所以实际操作中我们并不真的根据随机的q(x)确定每次对哪一维进行变换，而是通过一个固定的顺序进行模拟，如按1,2,3,…,n−1,n,1,2,…n,1,…这样的周期进行采样（这里需要强调的是千万不要把这里的周期跟马尔可夫链的周期混淆，这里的周期其实本质上是对q(x)这个均匀分布的模拟），即在时间k∗n+i(k≥0)，只有第i维上能够发生变化。通过这种近似得到的分布与原目标分布的近似度是非常高的，具体证明见[1]
现在，假如我们通过Gibbs抽样方法对Gibbs分布进行抽样，会发现mcmc的另一个重要的性质。先回顾一下Gibbs分布：
其中E(x)为能量函数，Z为归一化系数，即∑xe−E(x)T。通过mcmc采样法，我们可以巧妙的避开对归一化系数Z的计算，因为我们在mcmc的计算过程中并不会直接用到π(x)本身，而是用到了π(yi|x−i)这样的条件概率。因为
可以看到，等式右边上下的归一化系数可以被抵消掉，这样就成功地解决了分布本身难以计算（这里主要是因为归一化系数难以计算）这一问题（其实绝大多数时候，对于类马尔科夫随机场的模型，我们都很难显式地对Z进行计算）。这里就是我在前文所说的mcmc的另一个motivation了（憋到现在终于把它点出来了）。
五、总结
这两天为了弄清楚mcmc和玻尔兹曼机，看了不少资料，其中有些资料让我走了一些弯路，所以写下这篇博客，一方面希望能够让别人少走弯路，另一方面也希望借此进一步加深对mcmc的理解。关于走过的弯路，这里不得不说一下《神经网络与机器学习》这本书，这本书说实话本身写的就不是非常accessible，中文翻译版的更是有些不忍直视，另外书中还有一些错误，印象比较深的一个是中文第三版370页中返回次数的均值应为pi1−pi，而不是书中的(1−p−1i)。我个人非常不建议将这本书作为第一次学习某些神经网络的途径，如果非要看的话，最好看英文版的。
此外，这里顺便提一下著名的优化神器模拟退火算法，其实模拟退火算法的核心思想就是逐步降低系统温度T，对满足温度T下的吉布斯分布（能量函数即为目标函数）进行mcmc采样，直到T降到一定值时（理论上T无限趋向于0时可以保证找到全局最优，但是要采样等待时间的期望也是无限大，因为除了最优解以外的所有会随机到的样本的概率都是0）得到的采样结果作为近似最优解时的样本。
PS: 本文中在一些数学符号的使用上可能缺乏严谨性，在不影响理解的情况下我就不再多花时间做调整了。
来源: http://garygu.xin/wordpress/index.php/2017/09/16/mcmc/
references
[1] Brémaud P. Markov chains: Gibbs fields, Monte Carlo simulation, and queues[M]. Springer Science & Business Media, 2013.
[2]Andrieu C, De Freitas N, Doucet A, et al. An introduction to MCMC for machine learning[J]. Machine learning, 2003, 50(1-2): 5-43.
转载请注明：江海志の博客 » 马尔可夫链蒙特卡洛(MCMC)采样详解
