R语言Metropolis Hastings采样和贝叶斯泊松回归Poisson模型

在本文中,我想向你展示如何使用R的Metropolis采样从贝叶斯Poisson回归模型中采样。

Metropolis-Hastings抽样算法是一类马尔科夫链蒙特卡洛(MCMC)方法

由Kaizong Ye,Liao Bao撰写

Metropolis-Hastings算法

其主要思想是生成一个马尔科夫链使其平稳分布为目标分布。这种算法最常见的应用之一是在贝叶斯统计中从后验密度中取样,这也是本文的目标。

该算法规定对于一个给定的状态Xt,如何生成下一个状态  有一个候选点Y,它是从一个提议分布 ,中生成的,根据决策标准被接受,所以链条在时间t+1时移动到状态Y,即Xt+1=Y或被拒绝,所以链条在时间t+1时保持在状态Xt,即Xt+1=Xt。

Metropolis 采样

在Metropolis算法中,提议分布是对称的,也就是说,提议分布  满足 

,所以Metropolis采样器产生马尔科夫链的过程如下。

  • 选择一个提议分布. 在选择它之前,了解这个函数中的理想特征。
  • 从提议分布g中生成X0。
  • 重复进行,直到链收敛到一个平稳的分布。
  • 从 生成Y.
  • 从Uniform(0, 1)中生成U。
  • 如果 , 接受Y并设置Xt+1=Y,否则设置Xt+1=Xt。这意味着候选点Y被大概率地接受.
  • 递增t.

贝叶斯方法

正如我之前提到的,我们要从定义为泊松回归模型的贝叶斯中取样。

对于贝叶斯分析中的参数估计,我们需要找到感兴趣的模型的似然函数,在这种情况下,从泊松回归模型中找到。

现在我们必须为每个参数β0和β1指定一个先验分布。我们将对这两个参数使用无信息的正态分布,β0∼N(0,100)和β1∼N(0,100) 。

最后,我们将后验分布定义为先验分布和似然分布的乘积。

使用Metropolis采样器时,后验分布将是目标分布。


R语言中的block Gibbs吉布斯采样贝叶斯多元线性回归

阅读文章


计算方法

这里你将学习如何使用R语言的Metropolis采样器从参数β0和β1的后验分布中采样。

数据

首先,我们从上面介绍的泊松回归模型生成数据。


随时关注您喜欢的主题


似然函数

现在我们定义似然函数。在这种情况下,我们将使用这个函数的对数,这是强烈建议的,以避免在运行算法时出现数字问题。

先验分布

接下来我们定义参数β0和β1的先验分布。与似然函数一样,我们将使用先验分布的对数。

后验分布

由于我们是用对数工作的,我们把后验分布定义为似然函数的对数与先验分布的对数之和。

记住,这个函数是我们的目标函数f(.),我们要从中取样。

提议函数

最后,我们定义提议分布g(.|Xt)。由于我们将使用Metropolis采样器,提议分布必须是对称的,并且取决于链的当前状态,因此我们将使用正态分布,其平均值等于当前状态下的参数值。

Metropolis 采样器

最后,我们编写代码,帮助我们执行Metropolis采样器。在这种情况下,由于我们使用的是对数,我们必须将候选点Y被接受的概率定义为。

在这种情况下,我们将在算法的每20次迭代中为我们的最终链选择一个值。

由于MCMC链具有很强的自相关,它可能产生的样本在短期内无法代表真实的基础后验分布。那么,为了减少自相关,我们可以只使用链上的每一个n个值来稀释样本。

在这里,你可以看到ACF图,它给我们提供了任何序列与其滞后值的自相关值。在这种情况下,我们展示了初始MCMC链的ACF图和对两个参数的样本进行稀释后的最终链。从图中我们可以得出结论,所使用的程序实际上能够大大减少自相关。

结果

在这一节中,我们介绍了由Metropolis采样器产生的链以及它对参数β0和β1的分布。参数的真实值由红线表示。

与glm()的比较

现在我们必须将使用Metropolis采样得到的结果与glm()函数进行比较,glm()函数用于拟合广义linera模型。

下表列出了参数的实际值和使用Metropolis采样器得到的估计值的平均值。

结论

从结果来看,我们可以得出结论,使用Metropolis采样器和glm()函数得到的泊松回归模型的参数β0和β1的估计值非常相似,并且接近于参数的实际值。另外,必须认识到先验分布、建议分布和链的初始值的选择对结果有很大的影响,因此这种选择必须正确进行。


可下载资源

关于作者

Kaizong Ye拓端研究室(TRL)的研究员。

本文借鉴了作者最近为《R语言数据分析挖掘必知必会 》课堂做的准备。

​非常感谢您阅读本文,如需帮助请联系我们!

 
QQ在线咨询
售前咨询热线
15121130882
售后咨询热线
0571-63341498