• 承接《马尔可夫链蒙特卡罗法》,上篇中描述了马尔可夫链蒙特卡罗法的本质 —— 以马尔可夫链为概率模型的蒙特卡罗方法,在马尔可夫链上进行随机游走获得样本,构造时使得马尔可夫链的平稳分布就是我们抽样的目标分布,当马尔可夫链燃烧期度过后,可以认为获得的样本是来自于目标分布 p(x)p(x)p(x) 的,从而可以用遍历均值作为数学期望的估计。
  • Metropolis-Hastings算法是马尔可夫链蒙特卡罗法的一种具体实现,其算法流程类似于接受-拒绝采样法,任意选择初始样本 x0x_0x0​ 后,对于迭代变量取值从 111 到 nnn,记 xi−1=xx_{i-1}=xxi−1​=x,按照建议分布 q(x,x′)q(x,x')q(x,x′) 随机抽取候选状态 x′x'x′,而后分别从 U(0,1)U(0,1)U(0,1) 中随机选取参数 uuu 以及计算接受分布的概率 α(x,x′)\alpha(x,x')α(x,x′),当 u≤α(x,x′)u\leq\alpha(x,x')u≤α(x,x′) 时接受该候选状态,令 xt=x′x_t=x'xt​=x′,否则令 xt=x.x_t=x.xt​=x.
  • Gibbs抽样基于满条件概率分布,从一个初始样本 x(0)=[x1(0),x2(0),⋯,xk(0)]Tx^{(0)}=\big[x^{(0)}_1,x^{(0)}_2,\cdots,x^{(0)}_k\big]^Tx(0)=[x1(0)​,x2(0)​,⋯,xk(0)​]T 出发进行迭代,每次迭代得到服从联合分布的样本 x(i)=[x1(i),x2(i),⋯,xk(i)]Tx^{(i)}=\big[x^{(i)}_1,x^{(i)}_2,\cdots,x^{(i)}_k\big]^Tx(i)=[x1(i)​,x2(i)​,⋯,xk(i)​]T,最终得到样本序列 {x(0),x(1),⋯,x(m),x(m+1),⋯,x(n)}.\{x^{(0)},x^{(1)},\cdots,x^{(m)},x^{(m+1)},\cdots,x^{(n)}\}.{x(0),x(1),⋯,x(m),x(m+1),⋯,x(n)}. 样本的具体迭代方式是针对每一个随机变量分量按照满条件分布进行的,记第 i−1i-1i−1 次迭代得到的样本为 x(i−1)x^{(i-1)}x(i−1),那么第 iii 次迭代时,首先对 x1(i)x^{(i)}_1x1(i)​ 按照以下满条件概率分布随机抽样:p(x1(i)∣x2(i−1),x3(i−1),⋯,xk(i−1))p(x^{(i)}_1|x^{(i-1)}_2,x^{(i-1)}_3,\cdots,x^{(i-1)}_k)p(x1(i)​∣x2(i−1)​,x3(i−1)​,⋯,xk(i−1)​)对于 2,3,⋯,k−12,3,\cdots,k-12,3,⋯,k−1 分量,按照以下的满条件概率分布随机抽样:p(xj(i)∣x1(i),x2(i),⋯,xj−1(i),xj+1(i−1),⋯,xk(i−1)),j=2,⋯,k−1p(x^{(i)}_j|x^{(i)}_1,x^{(i)}_2,\cdots,x^{(i)}_{j-1},x^{(i-1)}_{j+1},\cdots,x^{(i-1)}_k)~,~j=2,\cdots,k-1p(xj(i)​∣x1(i)​,x2(i)​,⋯,xj−1(i)​,xj+1(i−1)​,⋯,xk(i−1)​) , j=2,⋯,k−1上述满条件概率分布是一种异步迭代的方式,对于 xk(i)x^{(i)}_kxk(i)​ 分量,按照以下的满条件概率分布随机抽样:p(xk(i)∣x1(i),x2(i),⋯,xk−1(i))p(x^{(i)}_k|x^{(i)}_1,x^{(i)}_2,\cdots,x^{(i)}_{k-1})p(xk(i)​∣x1(i)​,x2(i)​,⋯,xk−1(i)​)至此得到了整个第 iii 次迭代的样本:x(i)=[x1(i),x2(i),⋯,xk(i)]Tx^{(i)}=\big[x^{(i)}_1,x^{(i)}_2,\cdots,x^{(i)}_k\big]^Tx(i)=[x1(i)​,x2(i)​,⋯,xk(i)​]T下一步同样地以燃烧期后的样本集合遍历均值作为数学期望的估计。
  • Metropolis-Hastings算法中,通常我们需要对多元变量联合分布进行抽样,但大多数情况中对于多元变量联合分布的抽样是难以进行的,一种简化的方法是类似Gibbs抽样中对每一变量的条件分布依次进行抽样,得到等效于对整个多元变量一次抽样的结果,这种策略变体下的MH算法被称为单分量MH算法( Single Component Metropolis-Hastings ).

文章目录

  • 单分量MH算法.
  • Gibbs抽样.

单分量MH算法.

  • 为了得到容量为 nnn 的样本集合 {x(0),x(1),⋯,x(m),x(m+1),⋯,x(n)}\{x^{(0)},x^{(1)},\cdots,x^{(m)},x^{(m+1)},\cdots,x^{(n)}\}{x(0),x(1),⋯,x(m),x(m+1),⋯,x(n)},单分量MH算法以下面的步骤实现一次抽样:记第 i−1i-1i−1 次迭代时第 jjj 变量样本值为 xj(i−1)x^{(i-1)}_jxj(i−1)​,在第 iii 次迭代时,我们按照建议分布 q(xj(i−1),xj∣x−j(i))q(x^{(i-1)}_j,x_j|x^{(i)}_{-j})q(xj(i−1)​,xj​∣x−j(i)​) 产生分量 xjx_jxj​ 的候选值 xj′(i)x'^{(i)}_jxj′(i)​,其中:x−j(i)=[x1(i),x2(i),⋯,xj−1(i),xj+1(i−1),⋯,xk(i−1)]Tx^{(i)}_{-j}=\Big[x^{(i)}_{1},x^{(i)}_{2},\cdots,x^{(i)}_{j-1},x^{(i-1)}_{j+1},\cdots,x^{(i-1)}_{k}\Big]^Tx−j(i)​=[x1(i)​,x2(i)​,⋯,xj−1(i)​,xj+1(i−1)​,⋯,xk(i−1)​]T
  • 下一步计算接受概率:α(xj(i−1),xj′(i)∣x−j(i))=min⁡{1,p(xj′(i)∣x−j(i))⋅q(xj′(i),xj(i−1)∣x−j(i))p(xj(i)∣x−j(i))⋅q(xj(i−1),xj′(i)∣x−j(i))}\alpha(x^{(i-1)}_j,x'^{(i)}_j|x^{(i)}_{-j})=\min\Big\{1,\frac{p(x'^{(i)}_j|x^{(i)}_{-j})\cdot q(x'^{(i)}_j,x^{(i-1)}_j|x^{(i)}_{-j})}{p(x^{(i)}_j|x^{(i)}_{-j})\cdot q(x^{(i-1)}_j,x'^{(i)}_j|x^{(i)}_{-j})}\Big\}α(xj(i−1)​,xj′(i)​∣x−j(i)​)=min{1,p(xj(i)​∣x−j(i)​)⋅q(xj(i−1)​,xj′(i)​∣x−j(i)​)p(xj′(i)​∣x−j(i)​)⋅q(xj′(i)​,xj(i−1)​∣x−j(i)​)​}并按照和一般MH算法相同的步骤来决定是否接受候选值 xj′(i)x'^{(i)}_jxj′(i)​,如果 U(0,1)U(0,1)U(0,1) 抽样的结果满足 u≤α(xj(i−1),xj′(i)∣x−j(i))u\leq\alpha(x^{(i-1)}_j,x'^{(i)}_j|x^{(i)}_{-j})u≤α(xj(i−1)​,xj′(i)​∣x−j(i)​) 则令 xj(i)=xj′(i)x^{(i)}_j=x'^{(i)}_jxj(i)​=xj′(i)​,否则令 xj(i)=xj(i−1).x^{(i)}_j=x^{(i-1)}_j.xj(i)​=xj(i−1)​. 其余分量在这个过程中都视为 xjx_jxj​ 的参数保持不变,对应的马尔可夫链的转移概率可以表示如下:p(xj(i−1),xj′(i)∣x−j(i))=α(xj(i−1),xj′(i)∣x−j(i))⋅q(x(i−1),xj′(i)∣x−j(i))p(x^{(i-1)}_j,x'^{(i)}_j|x^{(i)}_{-j})=\alpha(x^{(i-1)}_j,x'^{(i)}_j|x^{(i)}_{-j})\cdot q(x^{(i-1)},x'^{(i)}_j|x^{(i)}_{-j})p(xj(i−1)​,xj′(i)​∣x−j(i)​)=α(xj(i−1)​,xj′(i)​∣x−j(i)​)⋅q(x(i−1),xj′(i)​∣x−j(i)​)

  • 出于可视化直观性,下图是对于二元随机变量进行SCMH算法的迭代过程,熟悉数值优化算法的话可能会类似的图示感到熟悉 —— 坐标上升法的等高线图也呈现这样平行于坐标轴的迭代过程。因为在更新 x1,x2x_1,x_2x1​,x2​ 中某一个时,另一个都被视作参数而固定,最终的结果表现为只会在水平或垂直方向有移动,从而产生新的样本点,另外由于SCMH算法有 1−α(xj(i−1),xj′(i)∣x−j(i))1-\alpha(x^{(i-1)}_j,x'^{(i)}_j|x^{(i)}_{-j})1−α(xj(i−1)​,xj′(i)​∣x−j(i)​) 的概率拒绝候选值 xj′(i)x'^{(i)}_jxj′(i)​,因此有可能在某些相邻的时刻样本点不发生移动。

Gibbs抽样.

  • 我们考虑在单分量Metropolis-Hastings算法中定义建议分布为当前分量随机变量的满条件概率分布,如下:q(x,x′)=p(xj′∣x−j)q(x,x')=p(x'_j|x_{-j})q(x,x′)=p(xj′​∣x−j​)此时进行验证可以发现:α(x,x′)=min⁡{1,p(x′)⋅q(x′,x)p(x)⋅q(x,x′)}=min⁡{1,p(x′)⋅p(xj∣x−j′)p(x)⋅p(xj′∣x−j)}=min⁡{1,p(xj′∣x−j′)⋅p(x−j′)⋅p(xj∣x−j′)p(xj∣x−j)⋅p(x−j)⋅p(xj′∣x−j)}\begin{aligned}\alpha(x,x')&=\min\Big\{1,\frac{p(x')\cdot q(x',x)}{p(x)\cdot q(x,x')}\Big\}\\ &=\min\Big\{1,\frac{p(x')\cdot p(x_j|x'_{-j})}{p(x)\cdot p(x'_j|x_{-j})}\Big\}\\ &=\min\Big\{1,\frac{p(x'_j|x'_{-j})\cdot p(x'_{-j})\cdot p(x_j|x'_{-j})}{p(x_j|x_{-j})\cdot p(x_{-j})\cdot p(x'_j|x_{-j})}\Big\} \end{aligned}α(x,x′)​=min{1,p(x)⋅q(x,x′)p(x′)⋅q(x′,x)​}=min{1,p(x)⋅p(xj′​∣x−j​)p(x′)⋅p(xj​∣x−j′​)​}=min{1,p(xj​∣x−j​)⋅p(x−j​)⋅p(xj′​∣x−j​)p(xj′​∣x−j′​)⋅p(x−j′​)⋅p(xj​∣x−j′​)​}​观察 x−jx_{-j}x−j​ 的定义可以发现,实际上 x−jx_{-j}x−j​ 与 x−j′x'_{-j}x−j′​ 的表达形式完全一致,所以我们有:p(x−j)=p(x−j′);p(⋅∣x−j)=p(⋅∣x−j′)p(x_{-j})=p(x'_{-j})~;~p(\cdot|x_{-j})=p(\cdot|x'_{-j})p(x−j​)=p(x−j′​) ; p(⋅∣x−j​)=p(⋅∣x−j′​)因此我们有:p(xj′∣x−j′)=p(xj′∣x−j),p(xj∣x−j′)=p(xj∣x−j)p(x'_j|x'_{-j})=p(x'_j|x_{-j})~,~p(x_j|x'_{-j})=p(x_j|x_{-j})p(xj′​∣x−j′​)=p(xj′​∣x−j​) , p(xj​∣x−j′​)=p(xj​∣x−j​)所以得到接受概率:α(x,x′)=min⁡{1,1}=1\alpha(x,x')=\min\{1,1\}=1α(x,x′)=min{1,1}=1也就是说在Gibbs抽样中所有的候选值 x′x'x′ 都不会被拒绝,我们限制满条件概率分布 p(xj′∣x−j)>0p(x'_j|x_{-j})>0p(xj′​∣x−j​)>0 即可满足对应马尔可夫链不可约,该马尔可夫链的转移核函数也就是满条件概率分布 p(xj′∣x−j)p(x'_j|x_{-j})p(xj′​∣x−j​),所以我们说吉布斯抽样是特殊的单变量Metropolis-Hastings算法,它按照以满条件概率分布作为建议分布进行随机抽样,此时所有的候选值都不会被拒绝( 接受概率分布恒为111 ).
  • 吉布斯抽样的限制也在于满条件概率分布,当满条件概率分布易于抽样时,吉布斯抽样既简洁又有效;当不满足该性质时,就需要定义合适的建议分布使用MH算法。
  • 对应于上面给出的等高线图,吉布斯抽样的迭代过程与之类似,不同之处在于吉布斯抽样由于不会发生拒绝,所以相邻迭代之间不存在样本不移动的情况。

【ML】单分量Metropolis-Hastings算法与Gibbs抽样相关推荐

  1. Metropolis–Hastings算法

    1蒙特卡洛方法 蒙特卡罗方法也称统计模拟方法,是一种以概率统计理论为指导的数值计算方法.蒙特卡洛方法的基本思想是,当所求解问题是某种随机事件出现的概率,或者是某个随机变量的期望值时,通过某种" ...

  2. MCMC中的Metropolis–Hastings算法与吉布斯采样

    Metropolis–Hastings算法是一种具体的MCMC方法,而吉布斯采样(Gibbs Sampling)是Metropolis–Hastings算法的一种特殊形式.二者在机器学习中具有重要作用 ...

  3. R语言实现MCMC中的Metropolis–Hastings算法与吉布斯采样

    创建测试数据 第一步,我们创建一些测试数据,用来拟合我们的模型.我们假设预测变量和因变量之间存在线性关系,所以我们用线性模型并添加一些噪音. trueA <- 5trueB <- 0tru ...

  4. mh采样算法推导_科学网—MCMC中的Metropolis Hastings抽样法 - 张金龙的博文

    Metropolis Hastings抽样法示例 jinlongzhang01@gmail.com Metropolis Hasting(下面简称MH)是蒙特卡罗马尔科夫链中一种重要的抽样方法.本文简 ...

  5. 也谈MCMC方法与Gibbs抽样

    个人博客传送门:点击打开链接 MCMC,即传说中的Markov Chain Mento Carlo方法.其主要用于统计推理中进行模拟抽样,尤其在贝叶斯推理中有着非常广泛的应用.如算法模型的后验参数估计 ...

  6. Metropolis Hastings MCMC when the proposal and target have differing support

    参考: Metropolis Hastings MCMC when the proposal and target have differing support Examples of errors ...

  7. 每天一道算法练习题--Day25 第一章 --算法专题 --- ----------蓄水池抽样

    蓄水池抽样 问题描述 算法描述 相关题目 总结 力扣中关于蓄水池抽样问题官方标签是 2 道,根据我的做题情况来看,可能有三四道.比重算是比较低的,大家可以根据自己的实际情况选择性掌握. 蓄水池抽样的算 ...

  8. Gibbs抽样方法详解

    Gibbs抽样方法的作用: 积分,期望或者联合分布很难计算,通常情况下当前面三个问题是NP问题时才需要Gibbs Sampling. 不然的话,直接计算就可以了嘛,既准确又快速,干嘛还要Gibbs S ...

  9. python实现Metropolis采样算法实例

    马尔可夫链蒙特卡罗(Markov Chain Monte Carlo ,MCMC )方法主要算法有:Metropolis-Hastings (MH)算法.Metropolis采样算法,以及吉布斯采样方 ...

最新文章

  1. 告诉你一种精简、优化代码的方式
  2. Oracle数据库设计规范
  3. 进程环境详解(四)---getenv、putenv和setenv函数详解
  4. java基础知识的一些细节问题
  5. [原创]RedisDesktopManager工具使用介绍
  6. 创意油墨飞溅效果的绿树矢量素材
  7. 【缺陷检测】基于matlab GUI形态学PCB电路板缺陷检测【含Matlab源码 821期】
  8. 刘汝佳第二章习题(前四)
  9. idea常用22种快捷键,脱离鼠标,便捷开发,赶紧收藏
  10. 【Get深一度】信号处理(三)——3db带宽
  11. 【经典】产品人面试中的一些软回答~~
  12. django.db.utils.DataError: (1366, “Incorrect string value: ‘\\xE5\\x85\\xAD\\xE5\\x8D\\x83‘ for colu
  13. 基于RiskPariyBlackLitterman的因子择时
  14. java发送邮件连接超时,Java邮件超时和连接超时处理
  15. 报错 The type类名 is already defined
  16. 桌面和文件管理器右键卡顿几秒的解决办法
  17. Unity资源加载发布到移动端iphone/ipad
  18. java软件工程师自我评价_Java开发工程师岗位自我评价范文
  19. 中国Azure新数据中心(区域)正式商用
  20. Youtube的个人视频门户:Portal:Youtube director

热门文章

  1. 如何设置亚马逊code促销活动?
  2. 曲线救国使用图片url
  3. 国家二级计算机mysql_全国计算机等级考试二级MySQL练习软件
  4. Marvell以太网交换芯片-88E6390x-简介
  5. 文本挖掘在网络舆情信息分析中的应用_笔记
  6. 【Python中的】列表生成式和字典生成式以及内置函数
  7. 【转】Android 自己收集的开源项目和文章集合(持续更新至2018.12.17)
  8. Ajax: 一个建立Web应用的新途径 [转]
  9. C语言/C++常见习题问答集锦之哆啦A梦
  10. mysql主从 主机宕机_MySQL主从宕机的解决方法