在这篇文章中,我将对多元线性回归使用block的Gibbs采样,得出block的Gibbs采样所需的条件后验分布。然后,对采样器进行编码,并使用模拟数据对其进行测试 。
贝叶斯模型
假设我们有一个样本量的主题。贝叶斯多元回归假设该向量是从多元正态分布中提取的 ,通过使用恒等矩阵,我们假设独立的观察结果。
到目前为止,这与多元正态回归相同。则将概率最大化可得出以下解 :
贝叶斯模型是通过指定为一个先验分布得到 。在此示例中,我将在以下情况下使用 先验值
block Gibbs
在对采样器进行编码之前,我们需要导出Gibbs采样器的 每个参数的后验条件分布。
条件后验取更多的线性代数。
这是一个非常漂亮和直观的结果。条件后验的协方差矩阵是协方差矩阵的估计,
还要注意,条件后验是一个多元分布。因此,在Gibbs采样器的每次迭代中,我们从后验绘制出一个完整的矢量 。
模拟
我模拟的 结果向量
。
运行 Gibbs采样器 会生成对真实系数和方差参数的估计。运行了500,000次迭代。周期为100,000次,10次迭代。
以下是MCMC链的图,其中真实值用红线表示。
# 计算后验摘要统计信息
post_dist %>%
group_by(para) %>%
summarise(median=median(draw),
lwr=quantile(.025),
upr=quantile(.975)) %>%
# 合并汇总统计信息
post_dist <- post_dist %>%
left_join(post_sum_stats, by='param')
# 绘制MCMC链
ggplot(post_dist,aes(x=iter,y=dra)) +
geom_line() +
geom_hline(aes(yintercept=true_vals))
这是修整后参数的后验分布:
01
02
03
04
ggplot(post_dist,aes(x=draw)) +
geom_histogram(aes(x=draw),bins=50) +
geom_vline(aes(xintercept = true_vals))
似乎能够获得这些参数的合理后验估计。为了确保贝叶斯估计器正常工作,我对1,000个模拟数据集重复了此过程。
这将产生1,000组后验均值和1,000组95%置信区间。平均而言,这1000个后验均值应以真实值为中心。平均而言,真实参数值应在95%的时间的置信区间内。
以下是这些评估的摘要。
“估计平均值”列是所有1,000个模拟中的平均后验平均值。偏差百分比均小于5%。对于所有参数,95%CI的覆盖率约为95%。
扩展
我们可以对该模型进行许多扩展。例如,可以使用除正态分布外的其他分布来拟合不同类型的结果。 例如,如果我们有二元数据,则可以将其建模为:
然后在上放一个先验分布。这个想法将贝叶斯线性回归推广到贝叶斯GLM。
在本文中概述的线性情况下,可以更灵活地对协方差矩阵建模。相反,假设协方差矩阵是对角线且具有单个公共方差。这是多元线性回归中的同方差假设。如果数据是分类的(例如,每个受试者有多个观察结果),我们可以使用反Wishart分布来建模整个协方差矩阵。