哈喽，大家好~
在机器学习实践中，我们经常需要回答如下问题：
-
某个特征在模型中是否“真正有效”？换言之，它的系数是否显著非零？ -
两个模型（如不同特征集、不同架构）孰优孰劣？ -
数据量不大、信噪比一般时，如何稳定评估证据，而不仅仅是依赖p值与显著性水平？
贝叶斯检验提供了一种以证据为核心的方式：直接比较数据在两种假设下的概率，用比值（贝叶斯因子）刻画数据对两种假设的相对支持程度。
本篇文章，聚焦于一个典型场景：线性回归中对某个特征系数是否为零的检验。
它对特征选择、可解释性分析以及小样本/中等样本场景下的稳健推断都很有价值。
贝叶斯检验核心理论
贝叶斯学派的假设检验从以下对象出发：
-
先验： -
似然： -
后验： -
边际似然（模型证据）：
当我们比较两个假设/模型 与 时，核心量是贝叶斯因子：
其对数形式（更数值稳定）为：
-
若 （或 ），数据支持 的程度大于支持 。 -
证据强度的经验尺度（Jeffreys尺度）通常将 分段解释，例如 为“弱证据”， 为“中等证据”， 为“强证据”等（仅作参考）。
在贝叶斯框架下，检验还可以结合先验下的模型概率，通过后验赔率进行理性决策：
若我们没有偏好，令 ，则后验赔率即为贝叶斯因子。
一个非常有用的技巧是 Savage–Dickey 密度比（SDDR），在嵌套模型（如 约束 的某个参数取固定值）且先验在该参数上可分解且相容（compatible）时：
其中 是在 下被固定的参数值。例如检验 时，
于是 。这个公式可以大幅简化计算，因为只需计算先验和后验在一点的密度值。
在线性回归中的贝叶斯检验
考虑带高斯噪声的线性回归：
我们对权重向量 采用零均值各向同性高斯先验：
其中 控制先验强度（先验“尺度”或正则强度）。
在已知噪声方差 的情况下，后验仍为高斯分布（共轭）：
后验协方差：
后验均值：
这是经典的线性高斯模型推导。MAP解等价于岭回归，其正则项 。
边际似然（证据）可通过线性高斯积分得到：
因此对任意模型 ，其对数证据为：
现在我们进行变量重要性检验：
-
: 第1个特征的系数 。 -
: ，即完整模型。
注意， 是 的嵌套模型。我们可以有两种方式计算 ：
方式A：证据比（模型证据的比值）
-
在 下， 。 -
在 下，将第1列从 中去掉或设为0，得到 ，则 。 -
则
计算时应基于 Cholesky 分解来保证数值稳定：若 ，则 ，且 可用回代求解。
方式B：Savage–Dickey 密度比（SDDR）
在 下，先验对 因子化，且 是 的点假设，因此
因为先验 ，其在0处的密度是：
后验 的边际分布 仍是正态：
于是后验在0处的密度为：
从而
这个公式简单直观：当后验均值 偏离0越明显（相对于后验方差 越小），即 越大，贝叶斯因子越倾向支持 ；当后验方差 相较先验方差 显著收缩时（信息增益），也倾向于支持 。
两种方式（证据比与SDDR）在本设定下是等价的，我们会在代码中验证两者的数值一致性（允许小数值误差）。
后验预测
在 上的后验预测分布为：
这为我们提供了预测不确定性的刻画，用于校准与模型诊断。
完整案例
我们构造一个带相关特征的虚拟数据集， 个样本、 个特征。
真实权重向量 中仅有第1个特征有效（非零），其余为零或近零。
使用贝叶斯线性回归（高斯先验），检验 vs 。
实现细节与训练流程
数据生成：
-
设置 ，生成带相关性的特征，协方差为 AR(1)： ，取 ，再对列做标准化。 -
设真实权重 （也可加一些较小扰动）。 -
设 ，生成高斯噪声。
先验与正则：
-
选 （可调），即先验方差 。 -
岭回归正则强度 。
贝叶斯推断与数值稳定：
-
后验参数：用 Cholesky 分解求解 与 （不直接求逆）。 -
证据与贝叶斯因子：用 的 Cholesky 分解求 log evidence。 -
SDDR：只需要 和 ，从 和 中读取。
模型训练与对照：用 sklearn 的 Ridge 拟合（alpha = ），比对 Ridge 的系数与后验均值 的一致性。
import numpy as np
import scipy.stats as st
from scipy.linalg import cho_factor, cho_solve, cholesky
from sklearn.linear_model import Ridge
import matplotlib.pyplot as plt
import seaborn as sns
np.random.seed(42)
# 1) 生成数据：相关特征 + 标准化
n = 200
p = 6
rho = 0.6
# 生成协方差矩阵（AR(1)）
Sigma_x = rho ** np.abs(np.subtract.outer(np.arange(p), np.arange(p)))
# 生成特征 X ~ N(0, Sigma_x)
X = np.random.multivariate_normal(mean=np.zeros(p), cov=Sigma_x, size=n)
# 对列做标准化（零均值、单位方差）
X = (X - X.mean(axis=0)) / X.std(axis=0, ddof=1)
# 真实权重（第1个特征有效）
w_true = np.array([2.0, 0.0, 0.0, 0.0, 0.0, 0.0])
sigma = 1.0# 噪声标准差已知（这里我们“设定”已知，实践中可估计）
y = X @ w_true + np.random.normal(0, sigma, size=n)
# 2) 先验设置
tau = 1.5
tau2 = tau**2
sigma2 = sigma**2
lambda_ridge = sigma2 / tau2 # 岭回归中的alpha
# 3) 后验计算：S = (X^T X / sigma^2 + I / tau^2)^{-1}, m = S X^T y / sigma^2
XT_X = X.T @ X
XT_y = X.T @ y
A = XT_X / sigma2 + np.eye(p) / tau2 # p x p
# Cholesky 分解求解 A^{-1} @ b，比直接逆更稳定
c_factor = cho_factor(A, lower=True, check_finite=False)
S = cho_solve(c_factor, np.eye(p), check_finite=False) # 协方差
m = cho_solve(c_factor, XT_y / sigma2, check_finite=False) # 均值（也是MAP解）
# 4) 模型证据（H1 & H0）与贝叶斯因子（方式A：证据比）
def log_marginal_evidence(X, y, sigma2, tau2):
# 计算 log p(y|X, sigma2, tau2) = -0.5[n log(2π) + log|Σ| + y^T Σ^{-1} y]
# Σ = sigma2 I + tau2 X X^T
n = X.shape[0]
Sigma = sigma2 * np.eye(n) + tau2 * (X @ X.T) # n x n
# Cholesky for Σ
L = cholesky(Sigma, lower=True, check_finite=False)
# log|Σ| = 2 sum log diag(L)
logdet = 2.0 * np.sum(np.log(np.diag(L)))
# solve Σ^{-1} y
# 先解 L z = y, 再解 L^T α = z
z = np.linalg.solve(L, y)
alpha = np.linalg.solve(L.T, z)
quad = y @ alpha
return-0.5 * (n * np.log(2 * np.pi) + logdet + quad)
# H1: 完整模型
log_ev_H1 = log_marginal_evidence(X, y, sigma2, tau2)
# H0: 去掉第1列（索引0）的嵌套模型
X0 = X[:, 1:] # p-1个特征
log_ev_H0 = log_marginal_evidence(X0, y, sigma2, tau2)
log_BF_10_evidence = log_ev_H1 - log_ev_H0
# 5) 贝叶斯因子（方式B：SDDR）
m1 = m[0]
S11 = S[0, 0]
log_prior_at_0 = -0.5 * np.log(2 * np.pi * tau2) # N(0,tau2)在0处的 log 密度
log_post_at_0 = -0.5 * np.log(2 * np.pi * S11) - 0.5 * (m1**2) / S11
# BF10 = prior(0) / post(0)
log_BF_10_SDDR = log_prior_at_0 - log_post_at_0
# 6) 对照：Ridge回归的系数是否与后验均值一致？
ridge = Ridge(alpha=lambda_ridge, fit_intercept=False) # 我们的 X 与 y 已均值中心化，无需截距
ridge.fit(X, y)
w_ridge = ridge.coef_
diff_norm = np.linalg.norm(w_ridge - m)
print("后验均值（m）:", m)
print("Ridge 系数:", w_ridge)
print("二者差的范数:", diff_norm)
print("log BF10（证据比）:", log_BF_10_evidence)
print("log BF10（SDDR）:", log_BF_10_SDDR)
# 7) 贝叶斯因子随样本量增长
ns = np.linspace(20, n, 10, dtype=int) # 10个点
logBFs = []
for n_sub in ns:
X_sub = X[:n_sub]
y_sub = y[:n_sub]
log_ev_H1_sub = log_marginal_evidence(X_sub, y_sub, sigma2, tau2)
log_ev_H0_sub = log_marginal_evidence(X_sub[:, 1:], y_sub, sigma2, tau2)
logBFs.append(log_ev_H1_sub - log_ev_H0_sub)
logBFs = np.array(logBFs)
# 8) 不同tau的敏感性分析（对数刻度）
taus = np.logspace(-1, 1, 20) # 0.1 到 10
logBFs_tau = []
for t in taus:
t2 = t**2
log_ev_H1_t = log_marginal_evidence(X, y, sigma2, t2)
log_ev_H0_t = log_marginal_evidence(X0, y, sigma2, t2)
logBFs_tau.append(log_ev_H1_t - log_ev_H0_t)
logBFs_tau = np.array(logBFs_tau)
# 9) 后验预测评估（用留出集/或随机抽样简化）
# 为了展示，我们用8:2划分
n_train = int(0.8 * n)
X_train, y_train = X[:n_train], y[:n_train]
X_test, y_test = X[n_train:], y[n_train:]
# 重新计算训练集的后验
XT_X_tr = X_train.T @ X_train
XT_y_tr = X_train.T @ y_train
A_tr = XT_X_tr / sigma2 + np.eye(p) / tau2
c_factor_tr = cho_factor(A_tr, lower=True, check_finite=False)
S_tr = cho_solve(c_factor_tr, np.eye(p), check_finite=False)
m_tr = cho_solve(c_factor_tr, XT_y_tr / sigma2, check_finite=False)
# 后验预测：mean = X_test m_tr, var = sigma^2 + diag(X_test S_tr X_test^T)
pred_mean = X_test @ m_tr
pred_var = sigma2 + np.sum(X_test @ S_tr * X_test, axis=1)
pred_std = np.sqrt(pred_var)
# 10) 可视化分析
plt.figure(figsize=(14, 10))
# 图1：w1 的先验与后验密度对比
plt.figure(figsize=(7, 6))
grid_w1 = np.linspace(-4, 4, 400)
prior_w1_pdf = st.norm.pdf(grid_w1, loc=0, scale=tau)
post_w1_pdf = st.norm.pdf(grid_w1, loc=m1, scale=np.sqrt(S11))
ci_prior = np.array([st.norm.ppf(0.025, 0, tau), st.norm.ppf(0.975, 0, tau)])
ci_post = np.array([st.norm.ppf(0.025, m1, np.sqrt(S11)), st.norm.ppf(0.975, m1, np.sqrt(S11))])
plt.plot(grid_w1, prior_w1_pdf, color='magenta', lw=2,
label=r'先验: $\mathcal{N}(0, \tau^2)$', alpha=0.9)
plt.plot(grid_w1, post_w1_pdf, color='deepskyblue', lw=3,
label=r'后验: $\mathcal{N}(m_1, S_{11})$', alpha=0.9)
plt.fill_between(grid_w1, 0, prior_w1_pdf,
where=(grid_w1>=ci_prior[0]) & (grid_w1<=ci_prior[1]),
color='pink', alpha=0.3)
plt.fill_between(grid_w1, 0, post_w1_pdf,
where=(grid_w1>=ci_post[0]) & (grid_w1<=ci_post[1]),
color='cyan', alpha=0.25)
plt.axvline(0, color='orange', ls='--', lw=2,
label=r'零值 ($H_0$)')
plt.title(r'图1：$w_1$ 的先验与后验密度对比', fontsize=13)
plt.xlabel(r'$w_1$')
plt.ylabel(r'密度')
plt.legend()
plt.show()
# 图2：log BF10 随样本数 n
plt.figure(figsize=(7, 6))
plt.plot(ns, logBFs, marker='o', color='limegreen', lw=3,
label=r'$\log \mathrm{BF}_{10}$ (证据比)', alpha=0.9)
plt.axhline(0.0, color='gray', ls='--', lw=1)
plt.axhline(np.log(3), color='gold', ls='--', lw=2,
label=r'Jeffreys 阈值 $\approx \ln(3)$')
plt.axhline(np.log(10), color='red', ls='--', lw=2,
label=r'Jeffreys 阈值 $\approx \ln(10)$')
plt.title(r'图2：$\log \mathrm{BF}_{10}$ 随样本数 $n$ 的增长', fontsize=13)
plt.xlabel(r'样本数 $n$')
plt.ylabel(r'$\log \mathrm{BF}_{10}$')
plt.legend()
plt.show()
# 图3：后验预测校准（预测均值 vs 真实）
plt.figure(figsize=(7, 6))
plt.scatter(y_test, pred_mean, s=45, c='magenta', alpha=0.8,
label=r'预测均值 $\mathbb{E}[y \mid x]$',
edgecolors='white', linewidth=0.7)
# 画 95% 预测区间
for i in range(len(y_test)):
y_true = y_test[i]
y_hat = pred_mean[i]
y_lo = y_hat - 1.96 * pred_std[i]
y_hi = y_hat + 1.96 * pred_std[i]
plt.plot([y_true, y_true], [y_lo, y_hi],
color='deepskyblue', alpha=0.6, lw=2)
miny = min(y_test.min(), pred_mean.min()) - 1.0
maxy = max(y_test.max(), pred_mean.max()) + 1.0
plt.plot([miny, maxy], [miny, maxy], color='limegreen',
ls='--', lw=2, label=r'理想线: $y = \hat{y}$')
plt.xlim(miny, maxy)
plt.ylim(miny, maxy)
plt.gca().set_aspect('equal', adjustable='box')
plt.title(r'图3：后验预测校准与不确定性', fontsize=13)
plt.xlabel(r'真实值 $y$')
plt.ylabel(r'预测均值 $\mathbb{E}[y\mid x]$')
plt.legend()
plt.show()
# 图4：log BF10 对先验尺度 τ 的敏感性
plt.figure(figsize=(7, 6))
plt.semilogx(taus, logBFs_tau, marker='s', color='orange', lw=3,
label=r'$\log \mathrm{BF}_{10}$ vs. $\tau$', alpha=0.9)
plt.axhline(0.0, color='gray', ls='--', lw=1)
plt.axvline(tau, color='red', ls=':', lw=2,
label=r'当前 $\tau$')
plt.title(r'图4：先验尺度 $\tau$ 的敏感性分析', fontsize=13)
plt.xlabel(r'先验尺度 $\tau$（对数刻度）')
plt.ylabel(r'$\log \mathrm{BF}_{10}$')
plt.legend()
plt.show()
w₁ 的先验与后验密度对比：
数据更新后，我们对 的信念如何变化？零点（ ）还合理吗？
log BF₁₀ 随样本数增长：
随着样本数增加，数据对 相对于 的证据如何积累？
后验预测校准图（预测均值 vs 真实）与不确定性区间：
在留出测试集上，模型是否“校准良好”？预测不确定性是否合理？
log BF₁₀ 对先验尺度 τ 的敏感性：
不同先验强度下，贝叶斯因子的稳定性如何？结论是否“先验鲁棒”？
与频率学派检验的对比
与p值的差异：
-
p值是“在零假设为真”的条件下观测到当前或更极端统计量的概率。它不直接比较两个模型的“证据”。 -
贝叶斯因子直接比较 与 ，因此体现“奥卡姆剃刀”：复杂模型的参数体积（通过行列式项）会被自动惩罚。 -
在小样本场景或需要模型选择时，BF更契合实际需求。
与岭回归的关系：
-
本文的先验与岭回归的 L2 正则具有等价关系，MAP 解等于岭回归解。 -
贝叶斯推断提供了不确定性（后验协方差）与模型证据，而不仅仅是点估计。
g-prior 与闭式BF：
-
Zellner’s g-prior： ，在嵌套模型下有基于 的闭式贝叶斯因子表达式，非常适合变量选择。但需要注意 的选择或对其置超先验（如hyper-g prior）。
未知 ：
-
若 未知，可对其设置共轭逆伽马先验（与 的正态先验构成正态-逆伽马共轭族）。此时后验边际分布中的许多量将以Student-t形式出现，SDDR有相应调整。实现略复杂，但同样有闭式表达。
广义线性模型（GLM）与分类任务：
-
对于逻辑回归等非线性模型，后验不再有闭式；可以使用拉普拉斯近似、变分推断或MCMC进行贝叶斯检验与BF估计。 -
在A/B测试（Bernoulli–Beta）中，SDDR与闭式BF也常见且实现简单（Beta–Binomial模型）。
总结
全文，我们从理论出发推导了后验、边际似然与贝叶斯因子，展示了两种计算方式（证据比与SDDR）及其等价性；也讨论了与频率学派检验的差异与互补性，并指出了进一步优化与注意事项。
这一套流程不仅适用于线性模型，也能启发你在更复杂的机器学习模型（广义线性模型、核方法、深度网络）中应用贝叶斯思想进行模型选择与不确定性评估。
最后
