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

原文链接：

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

---

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

贝叶斯模型

假设我们有一个样本量的主题。 贝叶斯多元回归假设该向量是从多元正态分布中提取的 ，通过使用恒等矩阵，我们假设独立的观察结果。正式地，

​

到目前为止，这与 环境中看到的多元正态回归相同。 则将概率最大化可得出以下解 ：

​

贝叶斯模型是通过指定为一个先验分布得到 。在此示例中，我将在以下情况下使用 先验值

​

block Gibbs

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

​

条件后验

​取更多的线性代数。

​

这是一个非常漂亮和直观的结果。 条件后验的协方差矩阵是协方差矩阵的频繁估计，

​

因此，在Gibbs采样器的每次迭代中，我们从后验画出一个完整的矢量 。

模拟

我模拟的 结果向量

运行 Gibbs采样器 会生成对真实系数和方差参数的估计。运行了500,000次迭代。修整周期为100,000次，修整了10次迭代。

以下是MCMC链的图，其中真实值用红线表示。

# calculate posterior summary statistics (stats not used in rest of code) post_sum_stats<-post_dist %>% group_by(param) %>% summarise(median=median(draw), lwr=quantile(draw,.025), upr=quantile(draw,.975)) %>% mutate(true_vals=c(tb,tphi)) # merge on summary statistics post_dist <- post_dist %>% left_join(post_sum_stats, by='param') # plot MCMC Chains ggplot(post_dist,aes(x=iter,y=draw)) + geom_line() + geom_hline(aes(yintercept=true_vals, col='red'), show.legend=FALSE)+ facet_grid(param ~ .,scale='free_y',switch = 'y') + theme_bw() + xlab('Gibbs Sample Iteration') + ylab('MCMC Chains') + ggtitle('Gibbs Sampler MCMC Chains by Parameter')

​

这是 修整后参数的后验分布：

ggplot(post_dist,aes(x=draw)) + geom_histogram(aes(x=draw),bins=50) + geom_vline(aes(xintercept = true_vals,col='red'), show.legend = FALSE) + facet_grid(. ~ param, scale='free_x',switch = 'y') + theme_bw() + xlab('Posterior Distributions') + ylab('Count') + ggtitle('Posterior Distributions of Parameters (true values in red)')

​

似乎我们能够获得这些参数的合理后验估计。 为了确保贝叶斯估计器正常工作，我对1,000个模拟数据集重复了此练习。

这将产生1,000组后验均值和1,000组95％可信区间。平均而言，这1000个后验均值应以事实为中心。平均而言，真实参数值应在95％的时间的可信区间内。

以下是这些评估的摘要。

​

“估计平均值”列是所有1,000个模拟中的平均后验平均值。非常好。偏差百分比均小于5％。对于所有参数，95％CI的覆盖率约为95％。

扩展

我们可以对该模型进行许多扩展。例如，可以使用除正态分布外的其他分布来适应不同类型的结果。 例如，如果我们有二元数据，则可以将其建模为：

​

这个想法将贝叶斯线性回归推广到贝叶斯GLM。

在本文中概述的线性情况下，可以更灵活地对协方差矩阵建模。相反，假设协方差矩阵是对角线且具有单个公共方差。这是多元线性回归中的同方差假设。如果数据是聚类的（例如，每个受试者有多个观察结果），我们可以使用反Wishart分布来建模整个协方差矩阵。