这篇回答节选自我在专栏《机器学习中的数学:概率统计》中的文章。在这篇文章中,我们通过全面细致的分析MCMC全过程的原理,帮助大家深入理解接受概率。
也欢迎关注我的知乎账号 @石溪 ,将持续发布机器学习数学基础及Python数据分析编程应用等方面的精彩内容。
这里我们要一举解决最核心的关键问题:对于任意给定的目标分布 ,我们如何找到以他为唯一平稳分布的马尔科夫链,并且基于马尔科夫链采样的方法,实现对其的近似采样。
找这么一个马尔科夫链,本质上就是要找到他的转移概率矩阵 ,那么首先先确立一个思考路径:有没有什么条件,使得只要我们的转移矩阵 满足了,就意味着目标分布 就是转移矩阵 对应的马尔科夫链的平稳分布呢?
还真有这么一个条件,这个条件就是马尔科夫链的细致平稳条件。
设有马尔科夫链 ,状态空间为 ,转移概率矩阵为 ,以及我们关注的状态分布 ,对于任意状态 ,任意时刻 都满足: ,则称状态分布 满足马尔科夫链的细致平衡条件,该状态分布 就是马尔科夫链 的平稳分布。
我们提前说一下,在这个过程中我们可以把马尔科夫链状态转移概率 按照需要简写成 或 ,表达的都是同一个意思。
看得出,细致平稳条件是一个充分条件,这个条件很强,只要满足了他,就相当于找到了我们要的马尔科夫链。我们来证明一下,很简单:
如果 ,
我们同时对两边的等式求和: ,显然也是成立的。
观察等式的右边, 的取值是独立于参数 的,可以提出来: ,同时依据马尔科夫状态转移矩阵的归一性,可知: ,因此归结起来就有:
到了这一步,其实已经证明结束了,不过可能还有同学没有明白,那我就最后画一张图来彻底点破:
从图中可以知道,对于每一个 的具体取值:
的计算过程实际上就是一次状态分布(行向量)与状态转移矩阵第 列(列向量)点乘的过程,得到的是一个数值,值恰好是状态分布(行向量)的第 个元素,那么如果让 取遍每一个索引值,那么最终得到的就是
, , ,...,
用一个式子综合起来就是:
这不正是平稳分布的定义式嘛,这正说明了,满足细致平稳条件的目标分布 就是以矩阵 为状态转移矩阵的马尔科夫链的平稳分布。
并且如果此时这个转移概率矩阵对应的马尔科夫链满足不可约、非周期和正常返,那么进一步说明该状态分布 是马尔科夫链的唯一平稳分布。当然如果找到了这个转移矩阵 ,一般而言都是满足这个条件的。
因此,一旦找到了状态转移矩阵,就确定了马尔科夫链。
可是问题来了,感觉还是很难找到这个矩阵啊,例如,我们随便找一个状态转移矩阵 ,一般是无法满足细致平稳条件的,即: ,注意这里我们把 写作 ,是一个意思,这是为了我们后面公式中的写法统一。
同时也把离散型和连续型的马尔科夫链都统一了起来。对于 ,离散型马尔科夫链时,我们称之为转移概率,连续型马尔科夫链我们称之为转移概率密度,或称转移核,本质上是一回事儿,可以类比pmf和pdf的区别和联系。
为了让这个等式能够相等,满足细致平稳条件,我们给他加点料,不是 不相等吗?下面我们让他两边相等,我们来介绍Metropolis-Hastings采样方法具体是怎么做的。
在这个采样方法中,针对每一个转移概率 ,我们再定义一个接受概率 ,使得每一个我们最终要构造的状态转移矩阵中的每一项 都是 和 相乘的结果:
即目标马尔科夫矩阵的状态转移矩阵中的每一项都满足:
在建议马尔科夫链中,从状态 转移到状态 的概率是 ,而在目标马尔科夫链中,从状态 到状态 的转移概率是 ,很显然, 表示一个接受概率,他满足 ,那么相比于建议马尔科夫链,在目标马尔科夫链中,从状态 到状态 的转移概率都减小了,但是同时我们知道,马尔科夫链的概率转移矩阵要满足 ,所有的 都减小了,那总有一个要增大?谁增大呢,显然是 增大了。
我们的目标也是借力打力,试想我们要找的目标是状态转移概率矩阵 ,他的唯一稳态是目标分布 ,我们参考上一节接受-拒绝采样方法的思想精髓:既然这个 矩阵难找,我们能不能也找一个所谓的“建议矩阵” ,这个 满足任意性,这样我们就能够利用一个明确的、便于随机游走的另一个马尔科夫链,基于转移概率的建议矩阵 中的每一个具体项 和接受概率 来共同决定最终从状态 到状态 的状态转移概率,以达到如下效果:
第一:最终的状态转移概率 就是我们目标马尔科夫链上状态 到状态 的状态转移概率?也就是说 ?
第二:通过上述方式得到的 满足细致平稳条件,使得目标分布 恰是矩阵 对应的唯一平稳分布?
换句话说,原本要找一个状态转移概率 ,满足细致平稳条件:
而现在则是在建议转移概率矩阵 可以任选的前提下,找到一个接受概率 ,替换掉 ,满足下面的等式成立:
对于任意给定的建议转移概率矩阵 ,这个接受概率 都一定存在吗?问题最终就卡在了这。然后幸运的是,这个答案是:一定存在。马尔科夫链蒙特卡洛方法中的Metropolis-Hastings采样中明确了接受概率 的表达式:
也就是说我们任选一个建议矩阵 ,那么由这个矩阵的每一项和对应的接受概率相乘:
,用这所有相乘得到的结果所对应构成的新矩阵,就是以目标分布 为唯一稳态分布的马尔科夫链的转移概率矩阵 。
这听起来有点神乎其神啊,是这么回事儿吗?我们怎么证明这一点?很简单,只要证明细致平稳条件 成立即可
那么在实际采样操作中,这个如何体现?
显然,分两种情况:
当 ,即 时,有: ,
当 ,即 时,有:
,
归结起来就是:
当 时,
当 时,
有了这个等式,我们再来推导 是否满足细致平稳条件:
到了这里,回顾一下我们的初心,我们的目标是证明细致平稳条件 ,而此时,我们已经推导出了:
那么接下来,我们只要说明:
,再化简一点,就是证明:
这个等式成立
我们把之前的等式关系摆在这里,方便我们推导:
当 时,
当 时,
由于 和 是任意的参数,把 写作 ,把 写作 就有:
当 时,
当 时,
对于 这个待证明的等式,我们还是分情况讨论:
当 时,我们则需要证明 ,此时对上了第一种情况,即 成立。
当 时,我们则需要证明 ,此时对上了第二种情况, 也是成立的。
那么,这就说明了任选一个建议马尔科夫链概率转移矩阵 ,配合接受概率 , 和目标分布 是满足细致平稳条件的。
那么接下来就是具体操作了,怎么理解,或者说怎么实现基于 矩阵的马尔科夫链随机游走同时再叠加一个接受概率 呢?
再说细致一点,我们就讨论从任意状态 向任意状态 进行状态转移的过程及状态转移概率,单纯基于 矩阵,状态 向每一个其他状态 的概率是 ,其中 ,如何叠加新的接受概率 ?
打个比方,比如有3个状态, , , , , , ,依照随机游走的过程,当以 的概率选择往状态$2$转移时,还要最后面临一道考验:在此基础上,此时以0.8的概率真正选择转移到状态2,以1-0.8=0.2的概率选择“半途而废”回到原状态1,这两个动作叠加起来,就是在我们 进行随机游走的过程。如图所示,我们在途中只看状态1到状态2和状态3转移的局部过程:
在目标马尔科夫链的转移概率矩阵 中:
同时验证: ,满足归一性。
那么在随机游走的过程中,怎么模拟这个接受-拒绝的过程呢?假设我们目前处于状态 ,那就当我们已经按照转移概率矩阵 ,依概率决定要转移到状态 时,增加最后一步:生成一个满足 之间均匀分布的随机数 ,如果 ,则决定随机游走的下一个状态是 ,否则下一个状态仍然是 。
还记得我们上一讲当中模拟概率的技巧吗?我们先生成一个 之间满足均匀分布的随机数 ,如果 ,则决定接受,否则就拒绝。
最后一个问题就是我们应该构建一个怎么样的建议转移概率矩阵 ?其实这个矩阵 选取的原则很简单:
第一是:必须要满足马尔科夫概率转移矩阵的归一性,即 。
第二是:获取 时刻的采样点 后,依概率选取 时刻的采样点 的过程简单易行,方便用计算机模拟。
下面我们来解决这个问题:或者说当前处在某个特定的采样点 ,我们如何决定下一个随机游走的采样点 ?我们可以看做,状态空间 中的采样值 和 分别对应了状态 和状态 ,那么从 时刻的采样点 转移到 时刻的采样点 的概率值就是存放在矩阵 位置上,以此类推构造出整个矩阵 。
其中一种比较好的定义方法是:
这是什么意思?他指的是从 时刻的某个指定采样点 转移到 时刻的指定采样点 的概率等于在均值为 ,方差为1的整体分布中取得 的概率值,感觉其实应该写成pdf的,因为正态分布是一个连续型的分布,但是这里因为是计算机采样,实际上我们已经把包括目标分布以及这里的正态分布都给离散化了,所以写成pmf更符合实际操作的需要。
那么这么构造的转移概率矩阵 行不行,好不好呢?
首先:行不行?
,也就是整个正态分布概率密度曲线上(这里严格说是离散化后的分布列)所有的取值点的概率求和,显然等于1,满足归一化的要求,因此可行性满足了。如下图所示:
这幅图就是离散化的以 为均值,1为方差的正态分布图,红色的竖线就是 的取值,而黑色的竖线就是示意的几个 的概率取值,我们将建议转移概率矩阵 的 ,也就是 定义为了这个离散化的以 为均值,1为方差的正态分布中 的取值概率,那么由于所有的 的取值概率之和为1,因此就有了: ,这样就彻底说明了上面的观点。
第二:好不好?
如果当前处于采样点 ,依概率 决定下一个采样点 ,难不难?太容易了,这不就是基于 的正态分布模型进行一次采样吗?这样采出的点就是依概率 的新采样点。利用程序在一个正态分布中采样,想想我们都做过多少遍了?
第三:意外之喜。
如下图所示,我们发现,依据正态分布曲线关于均值 对称,因此对于任意两个 和 ,均值为 ,方差为1的正态分布在 处的概率值和均值为 ,方差为1的正态分布在 处的概率值显然是相等的,即:
,即相当于 ,那么接受概率的表达式 就可以化简为 。
好的,最后我们来归纳一下马尔科夫链蒙特卡洛方法中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
本站所有内容均为互联网搜索引擎提供的公开搜索信息,本站不存储任何数据与内容,任何内容与数据均与本站无关,如有需要请联系相关搜索引擎包括但不限于百度,google,bing,sogou 等
© 2025 tinynews.org All Rights Reserved. 百科问答小站 版权所有