R语言代写Gibbs抽样的贝叶斯简单线性回归仿真分析

R语言代写Gibbs抽样的贝叶斯简单线性回归仿真分析

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

贝叶斯分析的许多介绍使用相对简单的教学实例。虽然这可以很好地介绍贝叶斯原理，但将这些原则扩展到回归并不是直截了当的。

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

贝叶斯模型

假设我们观察到的数据 $（y_i，x_i）$ 的 $我在 {1， dots，n }$ 。我们的模型 $ÿ$ 是

$y_i sim N（ beta_0 + beta_1 x_i， phi）$

有兴趣的是推断 $beta_0$ 和 $beta_1$ 。如果我们在方差项上放置系数的正常先验和在方差项之前的反伽马，那么这个数据的完整贝叶斯模型可以写成：

$y_i | beta_0， beta_1， phi sim N（ beta_0 + beta_1 x_i， phi）$
$beta_0 | mu_0， tau_0 sim N（ mu_0， tau_0）$
$beta_1 | mu_1， tau_1 sim N（ mu_1， tau_1）$
$phi | alpha， gamma sim IG（ alpha， gamma）$

假设超参数 $mu_0， tau_0， mu_1， tau_1， alpha$ 并且 $伽玛$ 已知，后验可以写成比例常数，

$p（ beta_0， beta_1， phi | vec y） propto Big [ prod_ {i = 1} ^ np（y_i | beta_0， beta_1， phi） Big] p（ beta_0 | mu_0， tau_0）p（ beta_1 | mu_1， tau_1）p（ phi | alpha， gamma）$

括号中的术语是数据或可能性的联合分布。其他术语包括参数的联合先验分布（因为我们隐含地假定先前独立，联合先前因子）。

这被认为是熟悉的表达方式：

$后路 propto可能性次先于$ 。

随附的R代码的第0部分从该模型生成用于指定“真实”参数的数据。我们稍后将使用此数据估计贝叶斯回归模型，以检查我们是否可以恢复这些真实参数。

Gibbs采样器

从这个后验分布中得出，我们可以使用吉布斯采样算法。吉布斯采样是一种迭代算法，它根据每个感兴趣参数的后验分布产生样本。它是通过以下列方式从每个参数的条件后验顺序绘制来实现的：

可以证明，在适当的老化期后，1000次抽取的剩余部分是从后部分布中抽出的。这些样本不是独立的。绘制的顺序是在后部空间中的随机游走，并且空间中的每个步骤取决于先前的位置。通常也会使用稀疏期（这里没有这样做）。减薄10意味着我们每10次抽奖。这个想法是，每一个平局可能依赖于以前的平局，但不能作为依赖于10日以前的平局。

条件后验分布

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

它有助于从完全非标准化的后验开始：

$p（ beta_0， beta_1， phi | vec y） propto phi ^ { - n / 2} e ^ { - frac {1} {2 phi} sum_ {i = 1} ^ {n }（y_i - （ beta_0 + beta_1x_i））^ 2} e ^ { - frac {1} {2 tau_0}（ beta_0 - mu_0）^ 2} e ^ { - frac {1} {2 tau_1}（ beta_1 - mu_1）^ 2} phi ^ { - （ alpha + 1）e ^ { - frac { gamma} { phi}}}$

为了找到参数的条件后验，我们简单地从不包括该参数的联合后验中删除所有项。例如，常数项 $beta_0$ 具有条件后验：

$p（ beta_0 | phi， beta_1， mu_0， tau_0， vec y） propto e ^ { - frac {1} {2 phi} sum_ {i = 1} ^ n（y_i - （ beta_0 + beta_1x_i））^ 2} e ^ { - frac {1} {2 tau_0}（ beta_0 - mu_0）^ 2}$

同样的，

$p（ beta_1 | beta_0， phi， mu_1， tau_1， vec y） propto e ^ { - frac {1} {2 phi} sum_ {i = 1} ^ {n}（y_i - （ beta_0 + beta_1x_i））^ 2} e ^ { - frac {1} {2 tau_1}（ beta_1 - mu_1）^ 2}$

条件后验可以被识别为另一个反伽马分布，具有一些代数操作。

$begin {aligned} p（ phi | beta_0， beta_1， alpha， gamma， vec y）＆ propto phi ^ { - n / 2} phi ^ { - （ alpha + 1）e ^ { - frac { gamma} { phi}}} e ^ { - frac {1} {2 phi} sum_ {i = 1} ^ {n}（y_i - （ beta_0 + beta_1x_i））^ 2} \＆= phi ^ { - （ alpha + frac {n} {2} + 1）} exp { - frac {1} { phi} Big [ frac {1} { 2 phi} sum_ {i = 1} ^ {n}（y_i - （ beta_0 + beta_1x_i））^ 2 + gamma Big]} \＆= IG（shape = alpha + frac {n } {2}， rate = frac {1} {2 phi} sum_ {i = 1} ^ {n}（y_i - （ beta_0 + beta_1x_i））^ 2 + gamma） end {aligned }$

的条件后验 $beta_0$ 和 $beta_1$ 不容易辨认。但是如果我们愿意使用网格方法，我们真的不需要通过任何代数。

网格方法

考虑 $beta_0$ 。网格方法是一种非常强力的方式（在我看来）从其条件后验分布中进行采样。这种条件分布只是一个函数 $beta_0$ 。因此，我们可以评估某些 $beta_0$ 值的密度。在R表示法中，这可以是grid = seq（-10,10，by = .001）。这个序列是点的“网格”。

然后在每个网格点评估的条件后验分布告诉我们该抽取的相对可能性。

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

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

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

例如，在函数rb0cond（）和rb1cond（）中，我实际上计算了导出的条件后验分布的对数。然后，我从所有评估的最大值中减去每个评估，然后从对数标度中取代，然后进行标准化。这倾向于处理这样的数字问题。

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

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

首先，它的计算成本很高。如我们所做的那样，通过代数并希望得到一个已知的后验分布可以在计算上更有效率 $披$ 。

其次，网格方法需要指定网格点的区域。如果条件后验在我们指定的[-10,10]的网格区间之外有明显的密度怎么办？在这种情况下，我们不会从条件后验得到准确的样本。记住这一点并尝试宽网格间隔非常重要。因此，我们需要聪明地处理数值问题，例如接近Inf的数字和R中的-Inf值。

仿真结果

现在我们可以从每个参数的条件后验中进行采样，我们可以实现Gibbs采样器。这在随附的R代码的第2部分中完成。它编码上面R中概述的相同算法。

结果很好。下图显示了1000个Gibbs样本的序列（删除了老化抽取并且未实施细化）。红线表示我们模拟数据的真实参数值。第四个图显示了截距和斜率项的联合后验，红线表示轮廓。

总而言之，我们首先导出了一个表达式，用于参数的联合分布。然后我们概述了用于从后部绘制样本的Gibbs算法。在这个过程中，我们认识到Gibbs方法依赖于每个参数的条件后验分布的顺序绘制。因为 $披$ ，这是一个容易识别，已知的分布。对于斜率和截距项，我们决定使用网格方法绕过代数。这样做的好处是我们支持很多代数。成本增加了计算复杂性，在选择适当的网格范围时会出现一些反复试验和数值问题。

如果您有任何疑问，请在下面发表评论。
相关阅读:
621
Java里的日期和时间学习
 [置顶] 宏扩展和参数扫描
 android 按字母搜索
 使用Eclipse EE开发web项目
 免解压版的Mysql的启动脚本，并且执行导入（windows）
高焕堂《android从程序员到架构师之路》 YY讲坛直面大师学习架构设计
 Android 计时与倒计时
 poj 1564 Sum It Up | zoj 1711 | hdu 1548 (dfs + 剪枝 or 判重)
字符型驱动程序的结构框架
原文地址：https://www.cnblogs.com/tecdat/p/10750592.html

R语言代写Gibbs抽样的贝叶斯简单线性回归仿真分析

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

贝叶斯模型

Gibbs采样器

条件后验分布

网格方法

仿真结果

如果您有任何疑问，请在下面发表评论。