基于Tidy方法在R中手动构建贝叶斯线性回归的先验、似然与后验

核心问题分析 #

贝叶斯后验的核心逻辑是：后验 ∝ 似然 × 先验。你手动计算出现偏差的原因，大概率是似然函数形式错误、先验与似然的结合逻辑不对，或是未正确关联线性预测与观测值、忽略了正态似然的残差方差项。

关键前提回顾 #

已知条件整理：

• 先验参数：截距μ₀=37.3，斜率β₀=-5.3；已通过MAD转Z分数得到先验对应的标准差sd_intercept、sd_slope
• 似然假设：响应变量服从正态分布，即 $y \sim N(\alpha + \beta x, \sigma_{\text{likelihood}})$，其中$\sigma_{\text{likelihood}}$是stan_glm拟合出的残差标准差
• 目标：通过dplyr手动计算后验，匹配stan_glm给出的后验截距43.6、斜率-5.6

修正后的dplyr代码实现 #

假设new_cars_1990s包含自变量x（如mpg）和因变量y（如hp），以下是完整计算流程：

library(dplyr)
library(stats)

# 1. 填入已知的先验参数与转换后的标准差
prior_params <- tibble(
  mu_intercept = 37.3,
  sd_intercept = 2.5, # 替换为你的MAD转换结果
  mu_slope = -5.3,
  sd_slope = 0.8,     # 替换为你的MAD转换结果
  sigma_likelihood = 3.2 # 从stan_glm结果中提取的残差标准差
)

# 2. 手动计算后验分布的参数均值（匹配stan_glm的后验估计）
posterior_estimates <- new_cars_1990s %>%
  # 生成覆盖后验范围的参数网格，步长越小精度越高
  crossing(
    intercept = seq(35, 50, by = 0.1),
    slope = seq(-7, -4, by = 0.1)
  ) %>%
  # 计算每个参数组合的总似然（所有样本的似然乘积）
  group_by(intercept, slope) %>%
  mutate(
    linear_pred = intercept + slope * x,
    sample_likelihood = dnorm(y, mean = linear_pred, sd = prior_params$sigma_likelihood),
    total_likelihood = prod(sample_likelihood)
  ) %>%
  ungroup() %>%
  # 计算先验密度（截距与斜率的联合先验）
  mutate(
    prior_intercept = dnorm(intercept, mean = prior_params$mu_intercept, sd = prior_params$sd_intercept),
    prior_slope = dnorm(slope, mean = prior_params$mu_slope, sd = prior_params$sd_slope),
    total_prior = prior_intercept * prior_slope,
    # 后验密度 ∝ 似然 × 先验，无需归一化（加权均值仅需相对权重）
    posterior_density = total_likelihood * total_prior
  ) %>%
  # 计算后验均值（即stan_glm给出的点估计）
  summarize(
    posterior_intercept = weighted.mean(intercept, w = posterior_density),
    posterior_slope = weighted.mean(slope, w = posterior_density)
  )

# 查看结果
posterior_estimates

关键修正点说明 #

• 似然计算逻辑：必须基于每个样本的线性预测值intercept + slope * x，计算观测值y的正态密度，再求所有样本的似然乘积（样本独立假设）。之前的错误大概率是未正确构建线性预测或误用了似然形式。
• 先验与似然的结合：后验密度是总似然与联合先验的乘积，无需归一化——加权均值只需要权重的相对大小即可。
• 参数网格范围：seq的范围必须覆盖stan_glm给出的后验估计区间，确保参数网格包含后验分布的峰值区域。
• 残差标准差匹配：似然中的$\sigma_{\text{likelihood}}$必须与stan_glm拟合时使用的残差标准差完全一致，否则似然计算会出现系统性偏差。

验证与调整建议 #

如果结果仍有偏差：

1. 检查sigma_likelihood是否从stan_glm的summary()输出中正确提取（对应sigma项的估计值）
2. 确认MAD转SD的计算：正态分布中SD ≈ MAD / 0.6745，确保转换系数正确
3. 缩小参数网格的步长（如by = 0.01），提升计算精度

内容的提问来源于stack exchange，提问作者hachiko