R语言马尔可夫MCMC中的Metropolis Hastings,MH算法抽样(采样)法可视化实例
原文链接:http://tecdat.cn/?p=26324
介绍
Metropolis Hastings 算法是一种非常简单的算法,用于从难以采样的分布中生成样本。
假设我们要从分布 π 中进行采样,我们将其称为“目标”分布。为简单起见,我们假设 π是实线上的一维分布,尽管它很容易扩展到一维以上(见下文)。
MH 算法通过模拟马尔可夫链来工作,其平稳分布为 π。这意味着,从长远来看,来自马尔可夫链的样本看起来像来自 π的样本。正如我们将看到的,该算法非常简单和灵活。
MH算法
转移核
要实现 MH 算法,用户必须提供一个“转移核”QQ。转移核只是一种在 给定当前位置(例如 xx)的情况下随机移动到空间中新位置(例如 y)的方式。也就是说,Q 是给定 x 在 y 上的分布,我们将其写成 Q(y|x)。在许多应用中,Q将是一个连续分布,在这种情况下 Q(y|x) 将是 y 上的密度,因此∫Q(y|x)dy=1(对于所有 x)。
例如,从当前位置 x 生成新位置 y 的一种非常简单的方法是向 x添加一个 N(0,1) 随机数。即设置y=x+N(0,1),或者转移y|x~N(x,1)。所以
文章图片
这种在当前位置x加上一些随机数得到y的核,在实际中经常使用,被称为“随机游走”核。
MH算法
使用转移核 Q 从目标分布 π 中采样的 MH 算法包括以下步骤:
- 初始化,X1=x1 。
- 对于 t=1,2,…
- 从 Q(y|xt)中采样 y。将 y 视为 xt+1 的“建议”值。
- 计算
-
文章图片
- AA通常被称为“接受概率”。
- 以概率 A“接受”提议的值,并设置 xt+1=y。否则设置 xt+1=xt。
- Metropolis 算法请注意,上面给出的示例随机游走建议 Q 对于所有 x,y 满足 Q(y|x)=Q(x|y) 任何满足这一点的建议都称为“对称”。当 Q 是对称时,MH 算法中 A 的公式 简化为:
文章图片
示例 为了帮助理解 MH 算法,我们现在做一个简单的例子:我们实现算法以从指数分布中采样:
文章图片
当然,以其他方式从指数分布中采样会容易得多;我们只是用它来说明算法。
请记住,π 被称为“目标”分布,因此我们调用函数来计算 π
target
:现在我们实现 MH 算法,使用上面提到的简单正态随机游走转移核 Q。
这是代码:
x = rep(0,10000)
x\[1\] = 3#初始化;我任意地将其设置为3
for(i in 2:10000){if(){
x\[i\] = proposed_x# 以最小(1,A)的概率接受移动。
} else {
x\[i\] = current_x# 否则就 "拒绝 "移动,并留在原地。
}
}
运行此代码后,我们可以绘制马尔可夫链 x 访问的位置(有时称为轨迹图)。
文章图片
请记住,我们设计此算法是为了从指数分布中采样。这意味着(只要我们运行算法足够长的时间!)x 的直方图应该看起来像一个指数分布。在这里我们检查一下:
hist(x)lines
文章图片
x 中的值的直方图确实提供了与指数分布的紧密拟合。
结束语 MH 算法的一个特别有用的特性是,即使 只知道π 是一个常数,它也可以实现:也就是说,对于一些已知的 f,π(x)=cf(x) , 但未知常数 c。这是因为该算法仅通过比率
文章图片
依赖于π 。
这个问题出现在贝叶斯应用中,其中后验分布与先验概率成正比,但比例常数通常是未知的。因此,MH 算法对于从后验分布进行采样以执行难以解析的贝叶斯计算特别有用。
文章图片
最受欢迎的见解
1.使用R语言进行METROPLIS-IN-GIBBS采样和MCMC运行
2.R语言中的Stan概率编程MCMC采样的贝叶斯模型
3.R语言实现MCMC中的Metropolis–Hastings算法与吉布斯采样
4.R语言BUGS JAGS贝叶斯分析 马尔科夫链蒙特卡洛方法(MCMC)采样
5.R语言中的block Gibbs吉布斯采样贝叶斯多元线性回归
6.R语言Gibbs抽样的贝叶斯简单线性回归仿真分析
7.R语言用Rcpp加速Metropolis-Hastings抽样估计贝叶斯逻辑回归模型的参数
8.R语言使用Metropolis- Hasting抽样算法进行逻辑回归
【R语言马尔可夫MCMC中的Metropolis Hastings,MH算法抽样(采样)法可视化实例】9.R语言中基于混合数据抽样(MIDAS)回归的HAR-RV模型预测GDP增长回归的HAR-RV模型预测GDP增长")
推荐阅读
- 使用Go语言写一个Http|使用Go语言写一个Http Server的实现
- C语言数据结构详细解析二叉树的操作
- JSF Facelets语言
- win10恢复右下角语言栏
- Linux下跨语言调用C++实践
- 升级Win10后更改显示语言无效如何处理
- 学习语言的20大Duolingo替代方案(最佳语言学习软件)
- Android Kotlin—— 语言结合
- 教你如何安装windows7 64中文语言包
- c#语言使用Unity粒子系统制作手雷爆炸