暂无图片
暂无图片
暂无图片
暂无图片
暂无图片

R语言马尔可夫MCMC中的METROPOLIS HASTINGS,MH算法抽样(采样)法可视化实例

拓端数据部落 2022-04-26
401

原文链接:http://tecdat.cn/?p=26324 

介绍

Metropolis Hastings 算法是一种非常简单的算法,用于从难以采样的分布中生成样本。

假设我们要从分布 π 中进行采样,我们将其称为“目标”分布。为简单起见,我们假设 π是实线上的一维分布,尽管它很容易扩展到一维以上(见下文)。

MH 算法通过模拟马尔可夫链来工作,其平稳分布为 π。这意味着,从长远来看,来自马尔可夫链的样本看起来像来自 π的样本。正如我们将看到的,该算法非常简单和灵活。

相关视频:R语言中RStan贝叶斯分析示例

MH算法

转移核

要实现 MH 算法,用户必须提供一个“转移核”Q。转移核只是一种在 给定当前位置(例如 x)的情况下随机移动到空间中新位置(例如 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 的“建议”值。

  • 计算

  • A通常被称为“接受概率”。

  • 以概率 A“接受”提议的值,并设置 xt+1=y。否则设置 xt+1=xt。

  • Metropolis 算法

    请注意,上面给出的示例随机游走建议 Q 对于所有 x,y 满足 Q(y|x)=Q(x|y) 任何满足这一点的建议都称为“对称”。当 Q 是对称时,MH 算法中 A 的公式 简化为:

该算法的这种特殊情况,具有 Q 对称,首先由 Metropolis 等人在 1953 年提出,因此它有时被称为“Metropolis 算法”。

示例

为了帮助理解 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 访问的位置(有时称为轨迹图)。



点击标题查阅往期内容


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


左右滑动查看更多


01

02

03

04


请记住,我们设计此算法是为了从指数分布中采样。这意味着(只要我们运行算法足够长的时间!)x 的直方图应该看起来像一个指数分布。在这里我们检查一下:

hist(x)

lines

x 中的值的直方图确实提供了与指数分布的紧密拟合。

结束语

MH 算法的一个特别有用的特性是,即使 只知道π 是一个常数,它也可以实现:也就是说,对于一些已知的 f,π(x)=cf(x) , 但未知常数 c。这是因为该算法仅通过比率依赖于π 。

这个问题出现在贝叶斯应用中,其中后验分布与先验概率成正比,但比例常数通常是未知的。因此,MH 算法对于从后验分布进行采样以执行难以解析的贝叶斯计算特别有用。




点击文末“阅读原文”

获取全文完整资料


本文选自《R语言马尔可夫MCMC中的METROPOLIS HASTINGS,MH算法抽样(采样)法可视化实例》。


点击标题查阅往期内容

python贝叶斯随机过程:马尔可夫链Markov-Chain,MC和Metropolis-Hastings,MH采样算法可视化
Python贝叶斯推断Metropolis-Hastings(M-H)MCMC采样算法的实现
Metropolis Hastings采样和贝叶斯泊松回归Poisson模型
Matlab用BUGS马尔可夫区制转换Markov switching随机波动率模型、序列蒙特卡罗SMC、M H采样分析时间序列
R语言RSTAN MCMC:NUTS采样算法用LASSO 构建贝叶斯线性回归模型分析职业声望数据
R语言BUGS序列蒙特卡罗SMC、马尔可夫转换随机波动率SV模型、粒子滤波、Metropolis Hasting采样时间序列分析
R语言Metropolis Hastings采样和贝叶斯泊松回归Poisson模型
R语言贝叶斯MCMC:用rstan建立线性回归模型分析汽车数据和可视化诊断
R语言贝叶斯MCMC:GLM逻辑回归、Rstan线性回归、Metropolis Hastings与Gibbs采样算法实例
R语言贝叶斯Poisson泊松-正态分布模型分析职业足球比赛进球数
R语言用Rcpp加速Metropolis-Hastings抽样估计贝叶斯逻辑回归模型的参数
R语言逻辑回归、Naive Bayes贝叶斯、决策树、随机森林算法预测心脏病
R语言中贝叶斯网络(BN)、动态贝叶斯网络、线性模型分析错颌畸形数据
R语言中的block Gibbs吉布斯采样贝叶斯多元线性回归
Python贝叶斯回归分析住房负担能力数据集
R语言实现贝叶斯分位数回归、lasso和自适应lasso贝叶斯分位数回归分析
Python用PyMC3实现贝叶斯线性回归模型
R语言用WinBUGS 软件对学术能力测验建立层次(分层)贝叶斯模型
R语言Gibbs抽样的贝叶斯简单线性回归仿真分析
R语言和STAN,JAGS:用RSTAN,RJAG建立贝叶斯多元线性回归预测选举数据
R语言基于copula的贝叶斯分层混合模型的诊断准确性研究
R语言贝叶斯线性回归和多元线性回归构建工资预测模型
R语言贝叶斯推断与MCMC:实现Metropolis-Hastings 采样算法示例
R语言stan进行基于贝叶斯推断的回归模型
R语言中RStan贝叶斯层次模型分析示例
R语言使用Metropolis-Hastings采样算法自适应贝叶斯估计与可视化
R语言随机搜索变量选择SSVS估计贝叶斯向量自回归(BVAR)模型
WinBUGS对多元随机波动率模型:贝叶斯估计与模型比较
R语言实现MCMC中的Metropolis–Hastings算法与吉布斯采样
R语言贝叶斯推断与MCMC:实现Metropolis-Hastings 采样算法示例
R语言使用Metropolis-Hastings采样算法自适应贝叶斯估计与可视化
视频:R语言中的Stan概率编程MCMC采样的贝叶斯模型
R语言MCMC:Metropolis-Hastings采样用于回归的贝叶斯估计


文章转载自拓端数据部落,如果涉嫌侵权,请发送邮件至:contact@modb.pro进行举报,并提供相关证据,一经查实,墨天轮将立刻删除相关内容。

评论