乐趣区

关于数据挖掘:Metropolis-Hastings采样和贝叶斯泊松回归Poisson模型附代码数据

全文下载链接: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 采样器产生马尔科夫链的过程如下。

  1. 抉择一个提议散布. 在抉择它之前,理解这个函数中的现实特色。
  2. 从提议散布 g 中生成 X0。
  3. 反复进行,直到链收敛到一个安稳的散布。
  • 从 生成 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 采样用于回归的贝叶斯预计

退出移动版