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

原文出处：拓端数据部落公众号

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

Metropolis-Hastings算法

Metropolis-Hastings抽样算法是一类马尔科夫链蒙特卡洛（MCMC）方法，其主要思想是生成一个马尔科夫链使其平稳分布为目标分布。这种算法最常见的应用之一是在贝叶斯统计中从后验密度中取样，这也是本文的目标。

该算法规定对于一个给定的状态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语言的Metropolis采样器从参数β0和β1的后验分布中采样。

数据

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

n <- 1000 #  样本大小
J <- 2 # 参数的数量
X <- runif(n,-2,2) # 生成自变量的值
beta <- runif(J,-2,2) #生成参数的值
y <- rpois(n, lambda = lambda) # 生成因变量的值

似然函数

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

LikelihoodFunction <- function(param){beta0 <- param[1] beta1 <- param[2] lambda <- exp(beta1*X + beta0)# 对数似然函数loglikelihoods <- sum(dpois(y, lambda = lambda, log=T)) return(loglikelihoods)
}

先验分布

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

beta0prior <- dnorm(beta0, 0, sqrt(100), log=TRUE)beta1prior <- dnorm(beta1, 0, sqrt(100), log=TRUE)return(beta0prior + beta1prior) #先验分布的对数

后验分布

由于我们是用对数工作的，我们把后验分布定义为似然函数的对数与先验分布的对数之和。记住，这个函数是我们的目标函数f(.)，我们要从中取样。

提议函数

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

Metropolis 采样器

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

        # 创建一个数组来保存链的值chain[1, ] <- startvalue # 定义链的起始值for (i in 1:iterations){# 从提议函数生成YY <- ProposalFunction(chain[i, ]) # 候选点被接受的概率PosteriorFunction(chain[i, ]))# 接受或拒绝Y的决策标准 if (runif(1) < probability) {chain[i+1, ] <- Y}else{ chain[i+1, ] <- chain[i, ]

由于MCMC链具有很强的自相关，它可能产生的样本在短期内无法代表真实的基础后验分布。那么，为了减少自相关，我们可以只使用链上的每一个n个值来稀释样本。在这种情况下，我们将在算法的每20次迭代中为我们的最终链选择一个值。

startvalue <- c(0, 0) # 定义链条的起始值
#每20次迭代选择最终链的值
for (i in 1:10000){if (i == 1){cfinal[i, ] <- chain[i*20,]} else {cfinal[i, ] <- chain[i*20,]# 删除链上的前5000个值
burnIn <- 5000

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

结果

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

与glm()的比较

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

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

##       True value Mean MCMC       glm
## beta0  1.0578047 1.0769213 1.0769789
## beta1  0.8113144 0.8007347 0.8009269

结论

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

最受欢迎的见解

1.R语言多元Logistic逻辑回归应用案例

2.面板平滑转移回归(PSTR)分析案例实现

3.matlab中的偏最小二乘回归（PLSR）和主成分回归（PCR）

4.R语言泊松Poisson回归模型分析案例

5.R语言回归中的Hosmer-Lemeshow拟合优度检验

6.r语言中对LASSO回归，Ridge岭回归和Elastic Net模型实现

7.在R语言中实现Logistic逻辑回归

8.python用线性回归预测股票价格

9.R语言如何在生存分析与Cox回归中计算IDI，NRI指标

拓端tecdat|R语言Metropolis Hastings采样和贝叶斯泊松回归Poisson模型相关推荐

拓端tecdat|R语言逻辑回归（Logistic回归）模型分类预测病人冠心病风险
最近我们被客户要求撰写关于冠心病风险的研究报告,包括一些图形和统计输出. 相关视频:R语言逻辑回归(Logistic回归)模型分类预测病人冠心病风险逻辑回归Logistic模型原理和R语言分类预测冠 ...
拓端tecdat|R语言用LOESS(局部加权回归)季节趋势分解（STL）进行时间序列异常检测
最近我们被客户要求撰写关于LOESS(局部加权回归)的研究报告,包括一些图形和统计输出. 这篇文章描述了一种对涉及季节性和趋势成分的时间序列的中点进行建模的方法.我们将对一种叫做STL的算法进行研究, ...
拓端tecdat|R语言向量误差修正模型 (VECMs)分析长期利率和通胀率影响关系
最近我们被客户要求撰写关于向量误差修正模型的研究报告,包括一些图形和统计输出. 向量自回归模型估计的先决条件之一是被分析的时间序列是平稳的.但是,经济理论认为,经济变量之间在水平上存在着均衡关系,可以 ...
拓端tecdat|R语言线性回归和时间序列分析北京房价影响因素可视化案例
最近我们被客户要求撰写关于北京房价影响因素的研究报告,包括一些图形和统计输出. 目的房价有关的数据可能反映了中国近年来的变化: 人们得到更多的资源(薪水),期望有更好的房子人口众多独生子女政策: ...
R语言泊松回归(poisson)模型案例:基于robust包的Breslow癫痫数据集
R语言泊松回归(poisson)模型案例:基于robust包的Breslow癫痫数据集目录 R语言泊松回归(poisson)模型案例:基于robust包的Breslow癫痫数据集 #数据加载
R语言使用ggplot2函数可视化需要构建泊松回归模型的计数目标变量的直方图分布并分析构建泊松回归模型的可行性
R语言使用ggplot2函数可视化需要构建泊松回归模型的计数目标变量的直方图分布并分析构建泊松回归模型的可行性目录
使用R语言的BNLearn包实现贝叶斯网络
转载自:http://f.dataguru.cn/thread-301701-1-1.html 1. 加载程序包导入数据 library(bnlearn) #CRAN中有,可以直接用install. ...
R语言垃圾邮件分类--朴素贝叶斯（机器学习）
邮件分类练习–朴素贝叶斯思路数据导入数据处理构建训练集和测试集词云展示数据降维训练模型模型测试提升模型一.数据导入文件目录为:C:\Users\kelanj\Documents\ ...
使用R语言进行Metroplis-in-Gibbs采样和MCMC运行分析
全文链接:http://tecdat.cn/?p=12200 对于许多模型,例如逻辑模型,没有共轭先验分布.因此,吉布斯采样不适用(点击文末"阅读原文"获取完整代码数据). 这篇文 ...
C语言metropolis方法,详解R语言MCMC:Metropolis-Hastings采样用于回归的贝叶斯估计
MCMC是从复杂概率模型中采样的通用技术. 蒙特卡洛马尔可夫链 Metropolis-Hastings算法问题如果需要计算有复杂后验pdf p(θ| y)的随机变量θ的函数f(θ)的平均值或期望 ...

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