全文下载链接：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采样器，提议散布必须是对称的，并且取决于链的以后状态，因而咱们将应用正态分布，其平均值等于以后状态下的参数值。
点击题目查阅往期内容
R语言中的block Gibbs吉布斯采样贝叶斯多元线性回归
左右滑动查看更多
01
02
03
04
Metropolis 采样器
最初，咱们编写代码，帮忙咱们执行Metropolis采样器。在这种状况下，因为咱们应用的是对数，咱们必须将候选点Y被承受的概率定义为。
# 创立一个数组来保留链的值 chain[1, ] <- startvalue # 定义链的起始值 for (i in 1:iterations){ # 从提议函数生成Y Y <- 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的估计值十分类似，并且靠近于参数的理论值。另外，必须意识到先验散布、倡议散布和链的初始值的抉择对后果有很大的影响，因而这种抉择必须正确进行。
本文摘选 《 R语言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采样用于回归的贝叶斯预计