Savage-Dickey密度比深解：bayestestR贝叶斯因子是怎么算出来的（附代码实现）

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在零假设点上取两个密度对先验分布做密度估计取 $x0$ 处的密度值 $\pi(0)$对后验分布做密度估计取 $x0$ 处的密度值 $\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) # ② 区间零假设 BFROPE 式检验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 值后可对照下表解读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., 2010Heck, 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），仅供参考