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

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

贝叶斯分析的许多介绍都使用了相对简单的教学实例（例如根据伯努利数据给出成功概率的推理）。虽然这可以很好地介绍贝叶斯原理，但是将这些原理扩展到回归并不是直接的。

这篇文章将概述这些原理如何扩展到简单的线性回归。我将导出感兴趣参数的后验条件分布，给出用于实现Gibbs采样器的R代码，并提出所谓的网格方法。

贝叶斯模型

假设我们观察数据

对于我们的模型是

有趣的是推断和。如果我们将正态先验放在系数上，将反伽玛先验放在方差项上，则此数据的完整贝叶斯模型可以写为：

假设超参数是已知的，后验可以被写入到一个比例常数，

括号中的项是是数据的联合分布或概率。其他项包括参数的联合先验分布。

R代码从指定的“真实”参数模型生成数据。我们稍后将用这个数据估计贝叶斯回归模型来检查是否可以恢复这些真实的参数。

tphi<-rinvgamma(1, shape=a, rate=g)
tb0<-rnorm(1, m0, sqrt(t0) )
tb1<-rnorm(1, m1, sqrt(t1) )
tphi; tb0; tb1;y<-rnorm(n, tb0 + tb1*x, sqrt(tphi))

吉布斯采样器

要从这种后验分布中得出，我们可以使用Gibbs采样算法。吉布斯采样是一种迭代算法，从每个感兴趣的参数的后验分布生成样本。它通过以下方式从每个参数的条件后验分布依次得出的：

可以看出，剩下的1,000个样本是从后验分布中抽取的。这些样本不是独立的。每个步骤都取决于先前的位置。

条件后验分布

要使用Gibbs，我们需要确定每个参数的条件后验。

为了找到参数的条件后验，我们简单地删除不包含参数后验的所有项。例如，常数项条件后验：

相似地，

条件后验可以被认为是另一个逆伽马分布。

条件后验不那么容易识别。但是如果使用网格方法，我们不需要通过代数方法。

考虑网格方法。网格方法是暴力方法，从其条件后验分布进行抽样。这个条件分布只是一个函数。所以我们可以评估密度值。在R中，可以是grid = seq（-10，10，by = .001）。这个序列是点的“网格”。

那么在每个网格点评估的条件后验分布告诉我们这个抽样的相对概率。

然后，我们可以使用R中的sample（）函数从这些网格中抽取，抽样概率与网格处的密度评估成比例。

for(i in 1:length(p) ){p[i]<- (-(1/(2*phi))*sum( (y - (grid[i]+b1*x))^2 ))  + ( -(1/(2*t0))*(grid[i] - m0)^2)}draw<-sample(grid, size = 1, prob = exp(1-p/max(p)))

这在R代码的第一部分的函数rb0cond（）和rb1cond（）中实现。

使用网格方法时遇到数值问题是很常见的。由于我们在评估网格中未标准化的后验，因此结果可能会变得相当大或很小。可能会在R中产生Inf和-Inf值。

例如，在函数rb0cond（）和rb1cond（）中，我实际上评估了条件后验分布的对数。然后，我进行归一化和对数化。

我们不需要使用网格方法来从后验条件抽样，因为它来自已知的分布。

请注意，这种网格方法有一些缺点。

首先，这在计算上是复杂的。通过代数，希望得到一个已知的后验分布，从而在计算上更有效率。

其次，网格方法需要指定网格点的区域。如果条件后验在我们指定的[-10,10]的网格区间之外具有显着的密度，在这种情况下，我们不会从条件后验得到准确的样本。并且需要广泛的网格区间进行实验。所以，我们需要灵活地处理数字问题，例如在R中接近Inf和-Inf值的数字。

仿真结果

现在我们可以从每个参数的条件后验进行采样，我们可以实现Gibbs采样器。

iter<-1000
burnin<-101
phi<-b0<-b1<-numeric(iter)
phi[1]<-b0[1]<-b1[1]<-6

结果很好。下图显示了1000个吉布斯（Gibbs）样本的序列。红线表示我们模拟数据的真实参数值。第四幅图显示了截距和斜率项的后验联合分布，红线表示等高线。

z <- kde2d(b0, b1, n=50)
plot(b0,b1, pch=19, cex=.4)
contour(z, drawlabels=FALSE, nlevels=10, col='red', add=TRUE)

总结一下，我们首先推导了一个表达式，用于参数的联合分布。然后我们概述了从后验抽取样本的Gibbs算法。在这个过程中，我们认识到Gibbs方法依赖于每个参数的条件后验分布，这是容易识别的已知的分布。对于斜率和截距项，我们决定用网格方法来规避代数方法。这样做的好处是我们可以绕开很多代数运算。代价是增加了计算复杂性。

参考文献

1.matlab使用贝叶斯优化的深度学习

2.matlab贝叶斯隐马尔可夫hmm模型实现

3.R语言Gibbs抽样的贝叶斯简单线性回归仿真

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

5.R语言中的Stan概率编程MCMC采样的贝叶斯模型

6.Python用PyMC3实现贝叶斯线性回归模型

7.R语言使用贝叶斯层次模型进行空间数据分析

8.R语言随机搜索变量选择SSVS估计贝叶斯向量自回归（BVAR）模型

9.matlab贝叶斯隐马尔可夫hmm模型实现

拓端tecdat|R语言Gibbs抽样的贝叶斯简单线性回归仿真分析相关推荐

拓端tecdat|R语言逻辑回归（Logistic回归）模型分类预测病人冠心病风险
最近我们被客户要求撰写关于冠心病风险的研究报告,包括一些图形和统计输出. 相关视频:R语言逻辑回归(Logistic回归)模型分类预测病人冠心病风险逻辑回归Logistic模型原理和R语言分类预测冠 ...
拓端tecdat|R语言用LOESS(局部加权回归)季节趋势分解（STL）进行时间序列异常检测
最近我们被客户要求撰写关于LOESS(局部加权回归)的研究报告,包括一些图形和统计输出. 这篇文章描述了一种对涉及季节性和趋势成分的时间序列的中点进行建模的方法.我们将对一种叫做STL的算法进行研究, ...
拓端tecdat|R语言向量误差修正模型 (VECMs)分析长期利率和通胀率影响关系
最近我们被客户要求撰写关于向量误差修正模型的研究报告,包括一些图形和统计输出. 向量自回归模型估计的先决条件之一是被分析的时间序列是平稳的.但是,经济理论认为,经济变量之间在水平上存在着均衡关系,可以 ...
拓端tecdat|R语言线性回归和时间序列分析北京房价影响因素可视化案例
最近我们被客户要求撰写关于北京房价影响因素的研究报告,包括一些图形和统计输出. 目的房价有关的数据可能反映了中国近年来的变化: 人们得到更多的资源(薪水),期望有更好的房子人口众多独生子女政策: ...
R语言中的block Gibbs吉布斯采样贝叶斯多元线性回归
在这篇文章中,我将对多元线性回归做同样的事情.我将得出阻塞的Gibbs采样器所需的条件后验分布.然后我将对采样器进行编码并使用模拟数据对其进行测试. 一个贝叶斯模型假设我们有一个样本大小的科目.我 ...
R语言基于copula的贝叶斯分层混合模型的诊断准确性研究
介绍在对诊断测试准确性的系统评价中,统计分析部分旨在估计测试的平均(跨研究)敏感性和特异性及其变异性以及其他测量.灵敏度和特异性之间往往存在负相关,这表明需要相关数据模型.由于用户,分析在统计上具有 ...
机械学习与R语言--Naive Bayes 朴素贝叶斯在R语言中的实现
为什么天气预报说70%概率下雨?为什么垃圾短信垃圾邮件被自动归类?这一切的基础算法便是朴素贝叶斯理论(算法有很多,这仅是其中之一). 1.由贝叶斯理论到朴素贝叶斯(naive bayes) 理论的基础 ...
拓端tecdat荣获掘金社区入驻新人奖
2021年7月,由掘金发起了"入驻成长礼"颁奖活动.本次活动邀请到知名开发者.服务机构代表等业界人士. 据了解,掘金社区"新入驻创作者礼"主要对已经积累了一定历 ...
R语言与抽样技术学习笔记（Jackknife）
R语言与抽样技术学习笔记(Randomize,Jackknife,bootstrap) Jackknife算法 Jackknife的想法在我很早的一篇博客<R语言与点估计学习笔记(刀切法与最小二 ...
拓端tecdat荣获2022年度51CTO博主之星
相信技术,传递价值,这是51CTO每一个技术创作者的动力与信念,2022 年度,拓端tecdat 作为新锐的数据分析咨询公司,在51CTO平台上,不断的输出优质的技术文章,分享前沿创新技术,输出最佳生 ...

拓端tecdat|R语言Gibbs抽样的贝叶斯简单线性回归仿真分析