在与贝叶斯相关的马尔可夫链蒙特卡洛方法中，为什么可以最大化后验概率？ 第1页

这篇回答节选自我在专栏《机器学习中的数学：概率统计》中的文章。在这篇文章中，我们通过全面细致的分析MCMC全过程的原理，帮助大家深入理解接受概率。

也欢迎关注我的知乎账号 @石溪 ，将持续发布机器学习数学基础及Python数据分析编程应用等方面的精彩内容。

1.问题的目标

这里我们要一举解决最核心的关键问题：对于任意给定的目标分布 ，我们如何找到以他为唯一平稳分布的马尔科夫链，并且基于马尔科夫链采样的方法，实现对其的近似采样。

找这么一个马尔科夫链，本质上就是要找到他的转移概率矩阵 ，那么首先先确立一个思考路径：有没有什么条件，使得只要我们的转移矩阵 满足了，就意味着目标分布 就是转移矩阵 对应的马尔科夫链的平稳分布呢？

还真有这么一个条件，这个条件就是马尔科夫链的细致平稳条件。

2.平稳分布判定：细致平稳条件

设有马尔科夫链 ，状态空间为 ，转移概率矩阵为 ，以及我们关注的状态分布 ，对于任意状态 ，任意时刻 都满足： ，则称状态分布 满足马尔科夫链的细致平衡条件，该状态分布 就是马尔科夫链 的平稳分布。

我们提前说一下，在这个过程中我们可以把马尔科夫链状态转移概率 按照需要简写成 或 ，表达的都是同一个意思。

看得出，细致平稳条件是一个充分条件，这个条件很强，只要满足了他，就相当于找到了我们要的马尔科夫链。我们来证明一下，很简单：

如果 ，

我们同时对两边的等式求和： ，显然也是成立的。

观察等式的右边， 的取值是独立于参数 的，可以提出来： ，同时依据马尔科夫状态转移矩阵的归一性，可知： ，因此归结起来就有：

到了这一步，其实已经证明结束了，不过可能还有同学没有明白，那我就最后画一张图来彻底点破：

从图中可以知道，对于每一个 的具体取值：

的计算过程实际上就是一次状态分布（行向量）与状态转移矩阵第 列（列向量）点乘的过程，得到的是一个数值，值恰好是状态分布（行向量）的第 个元素，那么如果让 取遍每一个索引值，那么最终得到的就是

， ， ，...，

用一个式子综合起来就是：

这不正是平稳分布的定义式嘛，这正说明了，满足细致平稳条件的目标分布 就是以矩阵 为状态转移矩阵的马尔科夫链的平稳分布。

并且如果此时这个转移概率矩阵对应的马尔科夫链满足不可约、非周期和正常返，那么进一步说明该状态分布 是马尔科夫链的唯一平稳分布。当然如果找到了这个转移矩阵 ，一般而言都是满足这个条件的。

因此，一旦找到了状态转移矩阵，就确定了马尔科夫链。

可是问题来了，感觉还是很难找到这个矩阵啊，例如，我们随便找一个状态转移矩阵 ，一般是无法满足细致平稳条件的，即： ，注意这里我们把 写作 ，是一个意思，这是为了我们后面公式中的写法统一。

同时也把离散型和连续型的马尔科夫链都统一了起来。对于 ，离散型马尔科夫链时，我们称之为转移概率，连续型马尔科夫链我们称之为转移概率密度，或称转移核，本质上是一回事儿，可以类比pmf和pdf的区别和联系。

为了让这个等式能够相等，满足细致平稳条件，我们给他加点料，不是 不相等吗？下面我们让他两边相等，我们来介绍Metropolis-Hastings采样方法具体是怎么做的。

3.Metropolis-Hastings采样方法

在这个采样方法中，针对每一个转移概率 ，我们再定义一个接受概率 ，使得每一个我们最终要构造的状态转移矩阵中的每一项 都是 和 相乘的结果：

即目标马尔科夫矩阵的状态转移矩阵中的每一项都满足：

在建议马尔科夫链中，从状态 转移到状态 的概率是 ，而在目标马尔科夫链中，从状态 到状态 的转移概率是 ，很显然， 表示一个接受概率，他满足 ，那么相比于建议马尔科夫链，在目标马尔科夫链中，从状态 到状态 的转移概率都减小了，但是同时我们知道，马尔科夫链的概率转移矩阵要满足 ，所有的 都减小了，那总有一个要增大？谁增大呢，显然是 增大了。

我们的目标也是借力打力，试想我们要找的目标是状态转移概率矩阵 ，他的唯一稳态是目标分布 ，我们参考上一节接受-拒绝采样方法的思想精髓：既然这个 矩阵难找，我们能不能也找一个所谓的“建议矩阵” ，这个 满足任意性，这样我们就能够利用一个明确的、便于随机游走的另一个马尔科夫链，基于转移概率的建议矩阵 中的每一个具体项 和接受概率 来共同决定最终从状态 到状态 的状态转移概率，以达到如下效果：

第一：最终的状态转移概率 就是我们目标马尔科夫链上状态 到状态 的状态转移概率？也就是说 ？

第二：通过上述方式得到的 满足细致平稳条件，使得目标分布 恰是矩阵 对应的唯一平稳分布？

换句话说，原本要找一个状态转移概率 ，满足细致平稳条件：

而现在则是在建议转移概率矩阵 可以任选的前提下，找到一个接受概率 ，替换掉 ，满足下面的等式成立：

对于任意给定的建议转移概率矩阵 ，这个接受概率 都一定存在吗？问题最终就卡在了这。然后幸运的是，这个答案是：一定存在。马尔科夫链蒙特卡洛方法中的Metropolis-Hastings采样中明确了接受概率 的表达式：

也就是说我们任选一个建议矩阵 ，那么由这个矩阵的每一项和对应的接受概率相乘：

，用这所有相乘得到的结果所对应构成的新矩阵，就是以目标分布 为唯一稳态分布的马尔科夫链的转移概率矩阵 。

4.M-H采样方法的严格证明

这听起来有点神乎其神啊，是这么回事儿吗？我们怎么证明这一点？很简单，只要证明细致平稳条件 成立即可

那么在实际采样操作中，这个如何体现？

显然，分两种情况：

当 ，即 时，有： ，

当 ，即 时，有：

，

归结起来就是：

当 时，

当 时，

有了这个等式，我们再来推导 是否满足细致平稳条件：

到了这里，回顾一下我们的初心，我们的目标是证明细致平稳条件 ，而此时，我们已经推导出了：

那么接下来，我们只要说明：

，再化简一点，就是证明：

这个等式成立

我们把之前的等式关系摆在这里，方便我们推导：

当 时，

当 时，

由于 和 是任意的参数，把 写作 ，把 写作 就有：

当 时，

当 时，

对于 这个待证明的等式，我们还是分情况讨论：

当 时，我们则需要证明 ，此时对上了第一种情况，即 成立。

当 时，我们则需要证明 ，此时对上了第二种情况， 也是成立的。

那么，这就说明了任选一个建议马尔科夫链概率转移矩阵 ，配合接受概率 ， 和目标分布 是满足细致平稳条件的。

5.如何理解随机游走叠加接受概率

那么接下来就是具体操作了，怎么理解，或者说怎么实现基于 矩阵的马尔科夫链随机游走同时再叠加一个接受概率 呢？

再说细致一点，我们就讨论从任意状态 向任意状态 进行状态转移的过程及状态转移概率，单纯基于 矩阵，状态 向每一个其他状态 的概率是 ，其中 ，如何叠加新的接受概率 ？

打个比方，比如有3个状态， ， ， ， ， ， ，依照随机游走的过程，当以 的概率选择往状态$2$转移时，还要最后面临一道考验：在此基础上，此时以0.8的概率真正选择转移到状态2，以1-0.8=0.2的概率选择“半途而废”回到原状态1，这两个动作叠加起来，就是在我们 进行随机游走的过程。如图所示，我们在途中只看状态1到状态2和状态3转移的局部过程：

在目标马尔科夫链的转移概率矩阵 中：

同时验证： ，满足归一性。

6.接受-拒绝过程的设计

那么在随机游走的过程中，怎么模拟这个接受-拒绝的过程呢？假设我们目前处于状态 ，那就当我们已经按照转移概率矩阵 ，依概率决定要转移到状态 时，增加最后一步：生成一个满足 之间均匀分布的随机数 ，如果 ，则决定随机游走的下一个状态是 ，否则下一个状态仍然是 。

还记得我们上一讲当中模拟概率的技巧吗？我们先生成一个 之间满足均匀分布的随机数 ，如果 ，则决定接受，否则就拒绝。

7.建议转移概率矩阵 的设计

最后一个问题就是我们应该构建一个怎么样的建议转移概率矩阵 ？其实这个矩阵 选取的原则很简单：

第一是：必须要满足马尔科夫概率转移矩阵的归一性，即 。

第二是：获取 时刻的采样点 后，依概率选取 时刻的采样点 的过程简单易行，方便用计算机模拟。

下面我们来解决这个问题：或者说当前处在某个特定的采样点 ，我们如何决定下一个随机游走的采样点 ?我们可以看做，状态空间 中的采样值 和 分别对应了状态 和状态 ，那么从 时刻的采样点 转移到 时刻的采样点 的概率值就是存放在矩阵 位置上，以此类推构造出整个矩阵 。

其中一种比较好的定义方法是：

这是什么意思？他指的是从 时刻的某个指定采样点 转移到 时刻的指定采样点 的概率等于在均值为 ，方差为1的整体分布中取得 的概率值，感觉其实应该写成pdf的，因为正态分布是一个连续型的分布，但是这里因为是计算机采样，实际上我们已经把包括目标分布以及这里的正态分布都给离散化了，所以写成pmf更符合实际操作的需要。

那么这么构造的转移概率矩阵 行不行，好不好呢？

首先：行不行？

，也就是整个正态分布概率密度曲线上（这里严格说是离散化后的分布列）所有的取值点的概率求和，显然等于1，满足归一化的要求，因此可行性满足了。如下图所示：

这幅图就是离散化的以 为均值，1为方差的正态分布图，红色的竖线就是 的取值，而黑色的竖线就是示意的几个 的概率取值，我们将建议转移概率矩阵 的 ，也就是 定义为了这个离散化的以 为均值，1为方差的正态分布中 的取值概率，那么由于所有的 的取值概率之和为1，因此就有了： ，这样就彻底说明了上面的观点。

第二：好不好？

如果当前处于采样点 ，依概率 决定下一个采样点 ，难不难？太容易了，这不就是基于 的正态分布模型进行一次采样吗？这样采出的点就是依概率 的新采样点。利用程序在一个正态分布中采样，想想我们都做过多少遍了？

第三：意外之喜。

如下图所示，我们发现，依据正态分布曲线关于均值 对称，因此对于任意两个 和 ，均值为 ，方差为1的正态分布在 处的概率值和均值为 ，方差为1的正态分布在 处的概率值显然是相等的，即：

，即相当于 ，那么接受概率的表达式 就可以化简为 。

8.步骤和演示

好的，最后我们来归纳一下马尔科夫链蒙特卡洛方法中metropolis-Hastings采样的步骤，并实际用python进行演示:

对于目标采样分布 ：

第一步：随机选定一个起始点 ，指定燃烧期 和稳定期 。

第二步：开始采样，每一轮采样都以上一轮的采样值 为均值，方差为1，生成一个正态分布，然后在这个正态分布中依概率随机选取一个值 。

第三步：在 的均匀分布中随机生成一个数 ，并指定接收概率 ，如果 ，则本轮新的采样值为 ，否则本轮新的采样值仍为上一轮的 。

重复第二步~第三步采样过程 次，结束后，保留后 次采样结果作为目标分布的近似采样。

我们还是对之前使用过的目标分布： 进行采样，燃烧期采样个数设定为 ，最终实际保留的有效采样点的个数为 。

代码片段：

       import random from scipy.stats import norm import matplotlib.pyplot as plt import numpy as np import seaborn seaborn.set()  # 目标采样分布pi def pi(x):     return (0.3 * np.exp(-(x - 0.3) ** 2) + 0.7 * np.exp(-(x - 2.) ** 2 / 0.3)) / 1.2113   m = 10000  # 燃烧期样本数 N = 100000  # 实际保留的有效样本数 sample = [0 for i in range(m + N)]  # 采样数组  sample[0] = 2  # 任选一个起始点，选择默认的0也可以，效果一样 # 基于接受概率，在建议马尔科夫链上随机游走采样 for t in range(1, m + N):     x = sample[t - 1]     x_star = norm.rvs(loc=x, scale=1, size=1)[0]  # 生成下一时刻随机游走的点x*     alpha = min(1, (pi(x_star) / pi(x)))  # 接受概率          u = random.uniform(0, 1)  # 生成满足0~1之间均匀分布的随机数     if u < alpha:  # 接受-拒绝的过程         sample[t] = x_star     else:         sample[t] = x  x = np.arange(-2, 4, 0.01) plt.plot(x, pi(x), color='r')  # 实际目标分布 plt.hist(sample[m:], bins=100, normed=True, color='b', edgecolor='k', alpha=0.6)  # 实际分布的近似采样 plt.show()     

运行结果：

此内容节选自我的专栏《机器学习中的数学：概率统计》，前三节免费试读，欢迎订阅：

当然还有《机器学习中的数学（全集）》系列专栏，欢迎大家阅读，配合食用，效果更佳~

有订阅的问题可咨询微信：zhangyumeng0422