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

最新推荐文章于 2023-11-20 17:11:35 发布

weixin_33893473

最新推荐文章于 2023-11-20 17:11:35 发布

阅读量774

点赞数

原文链接：http://www.cnblogs.com/tecdat/p/10751627.html

版权

在这篇文章中，我将对多元线性回归做同样的事情。我将得出阻塞的Gibbs采样器所需的条件后验分布。然后我将对采样器进行编码并使用模拟数据对其进行测试。

一个贝叶斯模型

假设我们有一个样本大小的 $ñ$ 科目。我们观察 $n \次1$ 结果向量 $ÿ$ 。贝叶斯多元回归假设该向量是从多元正态分布中得出的，其中均值向量是 $X \公测$ 和协方差矩阵 $\ phi我$ 。这里 $X$ 是观察到 $n \ times p$ 的协变量矩阵。注意，该矩阵的第一列是标识。参数矢量 $\公测$ 的 $p \次1$ ， $\披$ 是一种常见的方差参数， $一世$ 是 $n \ times n$ 单位矩阵。通过使用单位矩阵，我们假设独立观察。从形式上看，

$y \ sim N_n（X \ beta，\ phi I）$

到目前为止，这与频率设置中看到的多元正态回归相同。假设 $X$ 是完整列排名，最大化可能性产生解决方案：

$\ hat {\ beta} =（X'X）^ { - 1} X'y$

通过为 $\公测$ 和指定先验分布来获得贝叶斯模型 $\披$ 。对于这个例子，我将使用平坦的，不正确的先验 $\公测$ 和反伽马先验 $\披$ ：

$\ beta \ sim f（\ beta）\ propto 1$

$\ phi \ sim IG（\ alpha，\ gamma）$

我们假设超级参数是简单的。

联合后验分布

关节后验分布与...成正比

$p（\ beta，\ phi | y，X）\ propto p（y | \ beta，\ phi，X）p（\ beta | X）p（\ phi | X）$

我们可以这样写，因为我们假设事先独立。那是，

$p（\ beta，\ phi | X）= p（\ beta | X）p（\ phi | X）$ 。

替换分布，

$p（\ beta，\ phi | y，X）\ propto \ phi ^ { - n / 2} e ^ { - \ frac {1} {2 \ phi}（y - X \ beta）'（y - X \ beta）} \ phi ^ { - （\ alpha + 1）} e ^ { - \ frac {\ gamma} {\ phi}}$

block吉布斯采样器

在对采样器进行编码之前，我们需要导出Gibbs采样器的组件 - 每个参数的后验条件分布。

条件后验 $\披$ 是通过降低因子而 $\披$ 不是从关节后部和重新排列来找到的。这种情况很容易，因为没有什么可以放弃的：

$\ begin {aligned} p（\ phi | \ beta，\ alpha，\ gamma，y，X）＆\ propto \ phi ^ { - n / 2} e ^ { - \ frac {1} {2 \ phi}（ y - X \ beta）'（y - X \ beta）} \ phi ^ { - （\ alpha + 1）} e ^ { - \ frac {\ gamma} {\ phi}} \\＆\ propto \ phi ^ { - （\ alpha + \ frac {n} {2} + 1）} e ^ { - \ frac {1} {\ phi} \ Big [\ frac {1} {2}（y - X \ beta）' （y - X \ beta）+ \ gamma \ Big]} \\＆= IG（shape = \ alpha + \ frac {n} {2}，rate = \ frac {1} {2}（y - X \ beta ）'（y - X \ beta）+ \ gamma）\ end {aligned}$

条件后验 $\公测$ 需要更多的线性代数。

$\ begin {aligned} p（\ beta | \ phi，\ alpha，\ gamma，y，X）＆\ propto e ^ { - \ frac {1} {2 \ phi}（y - X \ beta）'（y - X \ beta）} \\＆\ propto e ^ { - \ frac {1} {2 \ phi}（y'y - 2 \ beta'X'y + \ beta'X'X \ beta）} \\ ＆\ propto e ^ { - \ frac {1} {2 \ phi}（\ beta'X'X \ beta - 2 \ beta'X'y）} \\＆\ propto e ^ { - \ frac {1} {2 \ phi}（\ beta'X'X \ beta - 2 \ beta'（X'X）（X'X）^ { - 1} X'y）} \\＆\ propto e ^ { - \ frac {1} {2 \ phi}（\ beta - （X'X）^ { - 1} X'y）'（X'X）（\ beta - （X'X）^ { - 1} X'y） } \\＆= N_p（（X'X）^ { - 1} X'y，\ phi（X'X）^ { - 1}）\ end {aligned}$

这是一个非常漂亮和直观的结果。因为我们在参数向量上使用平坦先验，所以参数向量的条件后验以最大似然估计为中心 $\ hat {\ beta} =（X'X）^ { - 1} X'y$ 。条件后验的协方差矩阵是协方差矩阵的频率估计， $Cov（\ hat \ beta）= \ phi（X'X）^ { - 1}$