Savage-Dickey密度比深解：bayestestR贝叶斯因子是怎么算出来的（附代码实现）
【免费下载链接】bayestestR :ghost: Utilities for analyzing Bayesian models and posterior distributions 
项目地址: https://gitcode.com/gh_mirrors/ba/bayestestR
你是否好奇 R 语言中 bayestestR 包输出的**贝叶斯因子（Bayes Factor）**背后究竟藏着什么？本文将用大白话拆解 Savage-Dickey 密度比原理，带你从源码层面看懂贝叶斯因子是怎么算出来的，并附上可直接运行的代码实现，帮助新手快速掌握这一贝叶斯统计核心工具 🔍
一、什么是贝叶斯因子？为什么需要 Savage-Dickey 密度比？
贝叶斯因子（BF） 是衡量"数据支持哪个假设"的比值：BF > 1 表示数据更支持备择假设，BF < 1 表示数据更支持零假设。它相当于 p 值的双向替代品——既能给出"证据支持零假设"的结论，也能量化支持的程度。
但严格计算贝叶斯因子需要对边际似然做高维积分，计算代价极高。幸运的是，Wagenmakers 等人（2010）提出的 Savage-Dickey 密度比提供了一个优雅的"捷径"：
只需比较先验分布和后验分布在零假设点上的"密度值"，二者相除就是贝叶斯因子！
换句话说：零假设点（通常是 0）在后验中的密度相对先验"缩水"了多少，数据就有多反对零假设 📉
)
图：面板 C 即 Savage-Dickey 密度比示意图——虚线为先验密度，黄色实线为后验密度，蓝点与红点分别是两者在零假设点（0）处的密度值，二者的比值就是贝叶斯因子。
二、Savage-Dickey 密度比原理：3 个关键步骤
步骤 1：用"密度"代替"概率"
零假设是"参数恰好等于 0"。连续分布中，单点的概率恒为 0，没法直接比较。于是我们改用概率密度：密度越高，说明该值在分布中越"可信"。
步骤 2：在零假设点上取两个密度
- 对先验分布做密度估计，取 $x=0$ 处的密度值 $\\pi(0)$
- 对后验分布做密度估计，取 $x=0$ 处的密度值 $\\pi(0\\mid D)$
步骤 3：两个密度相除
$$BF_{10} = \\frac{\\pi(0)}{\\pi(0\\mid D)}$$
- 若数据让后验"远离"了 0 → 后验在 0 处密度变低 → BF > 1（支持备择假设）
- 若数据让后验"聚焦"到 0 → 后验在 0 处密度升高 → BF < 1（支持零假设）
图：先验（虚线）与后验（黄色）密度对比。后验整体右移、变窄，其在零假设点处的密度相对先验大幅下降——这正是"数据反对零假设"的几何直观。
三、源码解析：bayestestR 是怎么一步步算出来的？
下面跟着源码走一遍完整流程（文件路径均相对项目根目录）。
1. 入口函数：智能分派
bayesfactor() 是统一入口，它会根据输入自动分派到合适的实现（见 R/bayesfactor.R 第 63–86 行）：
# R/bayesfactor.R（节选）
bayesfactor <- function(…, prior = NULL, direction = "two-sided",
null = 0, hypothesis = NULL, …) {
mods <- list(…)
if (length(mods) > 1) {
bayesfactor_models(…) # 多个模型 → 模型比较 BF
} else if (is.null(hypothesis)) {
bayesfactor_parameters( # 单参数 → Savage-Dickey 密度比
…, prior = prior, direction = direction, null = null
)
} else {
bayesfactor_restricted(…) # 带假设 → 有序限制 BF
}
}
2. 核心实现：logspline 密度估计 + 密度比
真正干活的是内部函数 .logbayesfactor_parameters()（R/bayesfactor_parameters.R 第 531–576 行）。简化后的核心逻辑如下：
# R/bayesfactor_parameters.R（核心逻辑，简化版）
.logbayesfactor_parameters <- function(posterior, prior, direction = 0, null = 0, …) {
if (length(null) == 1) {
# —— 点零假设：Savage-Dickey 密度比 ——
relative_loglikelihood <- function(samples) {
f_samples <- .logspline(samples) # ① logspline 拟合密度
d_samples <- logspline::dlogspline(null, f_samples, log = TRUE) # ② 取零假设点处的对数密度
if (direction < 0) {
norm_samples <- logspline::plogspline(null, f_samples) # ③ 单侧检验的归一化
} else if (direction > 0) {
norm_samples <- 1 – logspline::plogspline(null, f_samples)
} else {
norm_samples <- 1
}
d_samples – log(norm_samples)
}
} else {
# —— 区间零假设：比较区间内外的"相对可信度"变化 ——
# 用 logspline::plogspline() 计算区间概率，再做对数比值
}
# 关键一步：先验与后验的差值 = 对数贝叶斯因子
relative_loglikelihood(prior) – relative_loglikelihood(posterior)
}
逐行拆解 3 个要点：
| ① | .logspline(samples) | 用 logspline 算法（单调 B 样条）从 MCMC 样本中拟合平滑密度曲线，比简单核密度更精准 |
| ② | dlogspline(null, …) | 提取零假设点（默认 0）处的对数密度 |
| ③ | plogspline(…) | 单侧检验（direction = "left"/"right"）时，把非零一侧的概率归一化为 1，实现有序限制 |
最后第 575 行 relative_loglikelihood(prior) – relative_loglikelihood(posterior) 一行完成"先验密度比后验密度"——在对数空间相除变成相减，数值更稳定。函数返回的是 log_BF，用 as.numeric() 可转回普通 BF 值。
多参数模型则在外层 bayesfactor_parameters.data.frame()（第 475–484 行）中逐列循环调用上述核心函数，输出每个参数一行的 BF 表。
四、动手实践：3 行代码算出贝叶斯因子 💡
只需先验和后验的 MCMC 样本（向量、数据框或 stanreg/brmsfit 模型对象均可）：
library(bayestestR)
# 模拟 1000 个先验 / 后验样本
set.seed(123)
prior <- distribution_normal(1000, mean = 0, sd = 1)
posterior <- distribution_normal(1000, mean = 0.5, sd = 0.3)
# ① 点零假设 BF（默认 null = 0，即 Savage-Dickey 密度比）
bayesfactor_parameters(posterior, prior = prior, verbose = FALSE)
# ② 区间零假设 BF（ROPE 式检验，null 为区间）
bayesfactor_parameters(posterior, prior = prior,
null = c(-0.1, 0.1), verbose = FALSE)
# ③ 单侧（方向性）BF
bayesfactor_parameters(posterior, prior = prior,
direction = "right", verbose = FALSE)
bayesfactor_parameters()、bayesfactor_pointnull()、bayesfactor_rope() 是同一函数的三种默认配置（见 R/bayesfactor_parameters.R 第 195–241 行）。
图：主流工具（如 JASP）输出的贝叶斯 ANOVA 结果表，BF₁₀ 列的数值正是由上述 Savage-Dickey / 边际似然方法计算得到。
五、结果解读清单：贝叶斯因子数值怎么看？
输出结果默认是 log_BF（对数 BF），用 as.numeric() 提取 BF 值后，可对照下表解读：
| < 1/3 | 弱证据支持零假设 | 数据让零假设点密度升高 |
| 1/3 ~ 1 | 可忽略 | 证据极弱，两边都不明显 |
| 1 ~ 3 | 轶事级证据反对零假设 | 轻微偏向备择假设 |
| 3 ~ 10 | 中等证据 | 有一定支持 |
| 10 ~ 30 | 强证据 | 数据明显支持备择假设 |
| 30 ~ 100 | 很强证据 | 后验已大幅偏离零假设点 |
| > 100 | 极强证据 | 后验密度在 0 处几乎"塌陷" |
📌 小技巧：log 尺度上，log_BF ≈ 1 约相当于 BF ≈ 2.7（中等证据起点），log_BF ≈ 2.3 约相当于 BF ≈ 10（强证据起点），读数更方便。
六、避坑指南：4 个新手最容易踩的坑 ⚠️
样本量要够：源码中明确警告（R/bayesfactor_parameters.R 第 468–473 行），当样本少于 40,000 时会提示"Bayes factors might not be precise"。密度估计在分布"尾部/峰部"对样本量很敏感，建议 MCMC 至少跑 4 万条后验样本。
先验必须正确提供：BF 衡量的是"先验 → 后验"的信念变化。若 prior = NULL，函数会警告并默认把后验当前验，结果毫无意义。对 stanreg/brmsfit 模型可用 unupdate() 自动生成纯先验样本（见 R/unupdate.R）。
它只是"近似"：Savage-Dickey 密度比是局部替代假设（在零假设点附近取先验）下的贝叶斯因子近似（Wagenmakers et al., 2010；Heck, 2019 讨论了回归参数场景的注意点）。若先验很宽（如 t 分布重尾），BF 可能受先验形状影响较大。
有序限制 BF 只适用于事先假设：bayesfactor_restricted() 用于"参数 A > B"这类预设的方向性假设，不能事后挑选比较对象（见 R/bayesfactor_restricted.R 第 3–4 行的注释提醒）。
七、延伸阅读：项目文件导航 📚
想深入源码与文档，可从以下路径入手（相对项目根目录）：
- R/bayesfactor_parameters.R —— Savage-Dickey 密度比核心实现
- R/bayesfactor_restricted.R —— 有序限制 / 方向性假设的 BF 计算
- R/bayesfactor_models.R —— 多模型边际似然比较
- R/estimate_density.R —— 密度估计封装（支持 kernel / logspline 等方法）
- R/unupdate.R —— 从后验模型"还原"先验样本
- vignettes/bayes_factors.Rmd —— 官方贝叶斯因子完整教程（含推导）
- man/bayesfactor_parameters.Rd —— 函数参考文档
- tests/testthat/test-bayesfactor_parameters.R —— 相关单元测试
想从源码获取并本地构建，可执行：
git clone https://gitcode.com/gh_mirrors/ba/bayestestR
总结：bayestestR 的贝叶斯因子计算并不神秘——本质就是"logspline 拟合先验与后验密度 → 取零假设点处的密度 → 相除"三步。理解了 Savage-Dickey 密度比，你就掌握了从 R 的 MCMC 样本直接推断假设证据强度的钥匙 🗝️
【免费下载链接】bayestestR :ghost: Utilities for analyzing Bayesian models and posterior distributions 项目地址: https://gitcode.com/gh_mirrors/ba/bayestestR
创作声明：本文部分内容由AI辅助生成（AIGC），仅供参考
