MCMC、蒙特卡洛近似和Metropolis算法简介
MCMC 是Markov Chain Monte Carlo 的简称,但在传统模拟中有一个很重要的假设是样本是独立的(independent samples),这一点在贝叶斯统计尤其是高纬度的模型中很难做到。所以MCMC的目的就是运用蒙特卡洛模拟出一个马可链(Markov chain)。
如今,概率建模风靡一时,但是当我第一次了解它时,总有一件事情困扰我。 许多贝叶斯建模方法都需要计算积分,而我看到的任何工作示例似乎都使用高斯或伯努利分布,原因很简单如果您尝试使用比这更复杂的方法,它将成为分析的噩梦。 将贝叶斯模型限制在“表现良好”的分布的小子集中,可能会极大地阻碍你对问题建模的能力,所以我们必须找到克服这一限制的方法。
蒙特卡洛近似
如果我不想分析计算某个讨厌的积分怎么办?可以使用蒙特卡洛近似。
我们知道,我们可以通过使用目标分布的样本值计算期望通过使用目标分布的样本值计算样本均值。为什么重要?那么,期望是什么呢?
连续随机变量的期望。同样的过程也适用于离散的情况,只要改变求和的积分。
这种估计积分的方法由中心极限定理提供了一些很好的保证。首先,这是期望的无偏估计,其次,我们可以计算估计的方差。
使用蒙特卡罗样本计算积分是非常好的,但是我们如何从目标分布中抽取样本呢?绘制高斯或均匀样本很容易,但np.random会让你失望。
画样本最简单的方法是使用逆CDF方法但这依赖于获得逆CDF函数它通常没有一个很好的解析形式只对一维随机变量有意义。
Metropolis算法是许多马尔可夫链蒙特卡洛(MCMC)采样方法的组成部分之一。 当您可以访问的只是目标分布的pdf时,它使我们能够绘制样本。
MCMC方法需要注意的是我们不再采用独立样本所以我们不能保证估计的方差如何随着样本数量的增加而减少。如果样本是独立的,中心极限定理告诉我们估计值的方差将与样本数量(N)成反比地减少。对于MCMC,我们可以通过将样本数量从N调整为N_eff来对其贴现。N_eff(几乎)总是小于N,与链中样本的相关性有关。
Metropolis采样
Metropolis算法的步骤如下:
1.从目标分布域或先前分布的域中均匀采样起点。
2.在那时pdf。
3.根据一些状态转换函数,从当前位置走一步,为新样本提出建议。
4.计算新的pdf值。
5.计算新pdf /旧pdf的值。
6.如果比率大于1,请接受该步骤。
7.如果比率小于1:1.采样一个统一的随机变量。2.如果比率大于随机样本,请接受该步骤。
8.否则拒绝该建议,将当前职位作为示例添加并采取新的步骤。
请注意,在5–8中描述的过程等效于以概率为min(1,p(new)/ p(old))的伯努利概率接受样本,记住这个后面会用到……
Metropolis为何起作用?
对于任何MCMC方法,我们都希望确保一种称为详细平衡或可逆性的属性。
详细平衡
如果π满足,则π是马尔可夫链(1)的平稳分布。 我们可以通过对等式两边求和来证明这一点。 如果我们可以保证详细的平衡,那么我们也知道我们正在从马尔可夫链的固定分布中取样,我们将其作为目标分布。
详细平衡的思想是,由于每次转移的概率“质量”从状态i到状态j的转移与从状态j到状态i的转移是相同的,因此在链的每次转移之后,我们保持在平稳分布 。
现在,让我们展示一下Metropolis算法如何满足这一条件……
为了找到满足详细平衡的p(i,j),我们首先提出任意的“跳跃”概率q(i,j),然后仅通过接受概率为α(的“跳跃”来获得p(i,j)。 i,j)。 当“跳转”被拒绝时,状态保持j = i。 这种“接受”的想法并不是Metropolis算法独有的,它存在于MCMC的大多数变体中。
跳跃概率推导
这取决于α是有效概率分布。 因此,α的有效形式为:
Metropolis-Hastings 跳跃概率
如果跳跃概率是对称的,则可以简化为:
否则,它可以保留为完整形式,称为Metropolis-Hasting MCMC。
现在我们可以保证详细的平衡,我们可以让马尔可夫链式接管。 如果马尔可夫链是遍历的(所有状态都是不可约的),那么在某个时候,该链将到达平稳分布,并且我们能够从目标分布中获取样本。
你可能还注意到,由于alpha是π(j)/π(i)的函数。 这样的结果意味着目标分布无需标准化。 当将Metropolis用于难以计算证据项的贝叶斯后验估计时,这特别有用。
实现的注意事项
Metropolis算法的通用版本称为“随机行走Metropolis”,其中建议的状态为当前状态,再加上均值为零且协方差矩阵为σ²I的多元高斯。σ应选择为足够使得足够多的样本被拒绝大。 这是为了确保对目标分布进行良好的探索。
要注意的第二件事是老化的概念。 在马尔可夫链到达平稳分布之前采集的样本应删除,因为它们在链收敛之前不能代表目标分布。 确定要删除的样本数量有些困难,但通常,要删除的样本为10–20%(或有效的样本为10–100)。
Numpy代码实现
def metropolis(pi, dims, n_samples, burn_in=0.1, var=1):theta_ = np.random.randn(dims)*var samples = []while len(samples) < n_samples:theta = theta_ + np.random.randn(dims)*varratio = pi(theta)/pi(theta_)if np.random.rand(1) < ratio:sample = thetatheta_ = thetasamples.append(sample)else:sample = theta_samples.append(sample) samples = np.array(samples)return samples[int(samples*burn_in):,:]
我们可以看到这在两个高斯函数的和上的表现(注意这是一个非正态分布)。
from scipy.stats import multivariate_normal
def make_pdf(mean1, mean2, cov1, cov2): pdf1 = multivariate_normal(mean1, cov1)pdf2 = multivariate_normal(mean2, cov2) def pdf(x):return pdf1.pdf(x) * pdf2.pdf(x) return pdfpdf = make_pdf([3, 3],[-1, -1],np.array([[1,0.1],[0.1,1]], dtype=float),np.array([[1,0.1],[0.1,1]], dtype=float),
)
samples = metropolis(pdf, 2, 1_000, 0.1)
上面的gif显示了算法是如何遍历分布的,偶尔会在分布的两种不同模式之间跳转。注意,这也突出了metropolis算法的一个弱点,它处理相对较差的多模型分布。这是由于要探索一种新模式,步骤必须足够大,以便从一种模式希望到另一种模式。这要么需要一个大的步长,要么需要一个模态紧密分布。修改如哈密顿量MCMC可以帮助解决这一问题,但一般来说,这是大多数MCMC方法的一个问题。
论文地址:https://arxiv.org/pdf/1504.01896.pdf
作者:Alexander Bailey
deephub翻译组
MCMC、蒙特卡洛近似和Metropolis算法简介相关推荐
- 蒙特卡洛算法_MCMC、蒙特卡洛近似和Metropolis算法简介
MCMC 是Markov Chain Monte Carlo 的简称,但在传统模拟中有一个很重要的假设是样本是独立的(independent samples),这一点在贝叶斯统计尤其是高纬度的模型中很 ...
- 马尔可夫链蒙特卡洛(一):Metropolis算法
0. 引言 马尔可夫链蒙特卡洛(Markov Chain Monte Carlo,MCMC)是一类从复杂概率分布中进行采样的方法.Metropolis算法是严格意义上第一个MCMC算法,它于1953年 ...
- MCMC抽样---Metropolis算法
MCMC抽样-Metropolis算法 马尔科夫链蒙特卡洛抽样方法可追溯到1953年N.Metropolis等人在研究原子和分子的随机性运动问题时所引入的随机模拟方法.该方法被命名为Metropoli ...
- 如何实现马尔可夫链蒙特卡罗MCMC模型、Metropolis算法?
原文链接:http://tecdat.cn/?p=2687 什么是MCMC,什么时候使用它? MCMC只是一个从分布抽样的算法. 这只是众多算法之一.这个术语代表"马尔可夫链蒙特卡洛&quo ...
- 悼念!蒙特卡洛Metropolis算法贡献者之一Arianna Rosenbluth逝世
点击上方"3D视觉工坊",选择"星标" 干货第一时间送达 作者丨陈彩娴 编辑丨AI科技评论 AI科技评论消息,洛杉矶当地时间12月28日,Metropolis算 ...
- PRML第十一章读书笔记——Sampling Methods 拒绝采样/重要性采样/采样重要性重采样/数据增广IP算法、Metropolis算法/MH算法/吉布斯、切片采样、混合MC、估计配分函数
(终于把第十章读完了,这一章应该相对轻松.但这两天状态有待调整,所以没咋认真读) 目录 11.1 Basic Sampling Algorithms P526 标准概率分布 P528 拒绝采样 P53 ...
- 数据结构与算法:算法简介
数据结构与算法:算法简介 雪柯 大工生物信息 提笔为写给奋进之人 已关注 你说呢 . shenwei356 等 70 人赞同了该文章 引用自算法图解,作者[美] Aditya Bhargava 译袁国 ...
- 曲线的生成算法实现_PCGPlanet1-地形生成算法简介
比较常用的地形生成算法有三种: 四叉树算法,GeoMipmap算法,移动立方体算法 目前市面游戏采用的方案基本都是以这三种算法为基础实现的,下面依次进行介绍 四叉树算法 很经典的算法,在没有GPU的时 ...
- 爬山算法和模拟退火算法简介(转)
源:爬山算法和模拟退火算法简介 一. 爬山算法 ( Hill Climbing ) 介绍模拟退火前,先介绍爬山算法.爬山算法是一种简单的贪心搜索算法,该算法每次从当前解的临近解空间中选择一个最优解作为 ...
最新文章
- 拼多多成立技术顾问委员会,陆奇将领导相关工作
- 第五讲 类的封装和类的继承
- ML之LGBMRegressor(Competition):2018年全国大学生计算机技能应用大赛《住房月租金预测大数据赛》——设计思路以及核心代码—191017再次更新
- Java8新特性总结 - 4.方法引用
- android8.0学习(1)---Android Treble 概述
- canvas绘图数学知识总结
- 高通audio数据到Speaker播放流程
- python怎么用input输入列表_Python - 根据列表内容验证用户输入的最佳方法是什么?...
- Yii在控制层中引入模版进行渲染的几种方式。
- 制作RPG游戏的部分核心代码分析
- 软件测试之开发转测试相比于其他测试人员有哪些优势
- java里如何继承一个类_java如何继承类
- Travis CI(持续集成)
- Aras Innovator: AML包
- java程序假死_分析一个常见的java多线程通信问题(假死现象)
- 正片叠底(Multiply)和滤色(Screen)是两种基本的混合模式
- 所有学java的女生都进来看看
- vscode连接模拟器运行flutter项目
- [导入]转帖------牛逼顿
- 卡西欧350计算机度分秒转换,卡西欧FX-4500PA计算器怎样将如:12.58244度转换成度分秒啊...