在PyMC中实现基于专家信息先验的贝叶斯模型对比遇阻

• 我在对比多个带信息先验的正态分布贝叶斯模型时遇到技术瓶颈，需要帮助。每个模型的均值和标准差都基于领域知识设置了信息先验（informative priors）。
• 场景目标：对比不同专家的估计值，每个专家的估计对应一个先验分布，需将这些估计纳入模型对比流程。
• 尝试过的方法：用scipy和PyMC实现，但未找到参考资料；知道MCMC采样可行，也想过利用先验共轭性，但均无进展。
• 当前困境：用PyMC实现后，发现LOO值比贝叶斯因子（Bayes factor）更难解释，附上尝试的代码：

import pymc as pm
import numpy as np
import arviz as az

np.random.seed(123)
# 生成模拟数据：均值40，标准差10，样本量100
data = np.random.normal(loc=40, scale=10, size=100)

# 定义模型A：均值和标准差使用较宽的正态先验
with pm.Model() as model_a:
    # 均值先验：正态分布，均值40，标准差5
    mean = pm.Normal('mean', mu=40, sigma=5)
    # 标准差先验：正态分布，均值10，标准差5
    sigma = pm.Normal('sigma', mu=10, sigma=5)
    
    # 似然函数：正态分布，观测数据为data
    likelihood = pm.Normal('likelihood', mu=mean, sigma=sigma, observed=data)
    # 运行MCMC采样，获取后验样本
    trace_a = pm.sample(5000)

# 定义模型B：均值和标准差使用更窄的正态先验（对应更确定的专家估计）
with pm.Model() as model_b:
    # 均值先验：正态分布，均值40，标准差1
    mean = pm.Normal('mean', mu=40, sigma=1)
    # 标准差先验：正态分布，均值10，标准差1
    sigma = pm.Normal('sigma', mu=10, sigma=1) 
    
    # 似然函数：正态分布，观测数据为data
    likelihood = pm.Normal('likelihood', mu=mean, sigma=sigma, observed=data)
    # 运行MCMC采样，获取后验样本
    trace_b = pm.sample(5000)

# 手动计算BIC的函数
def calculate_bic(log_likelihood, num_parameters, num_data_points):
    bic = -2 * log_likelihood + num_parameters * np.log(num_data_points)
    return bic

# 计算每个模型的对数似然
with model_a:
    pm.compute_log_likelihood(trace_a)

with model_b:
    pm.compute_log_likelihood(trace_b)

# 计算LOO值并对比模型
a_loo = az.loo(trace_a)
df_comp_loo = az.compare({"A": trace_a, "B": trace_b})
az.plot_compare(df_comp_loo, insample_dev=False);

# 计算BIC
log_likelihood_a = model_a.logp().sum()
log_likelihood_b = model_b.logp().sum()

num_parameters_a = len(model_a.free_RVs)
num_parameters_b = len(model_b.free_RVs)

num_data_points = len(data)

bic_a = calculate_bic(log_likelihood_a, num_parameters_a, num_data_points)
bic_b = calculate_bic(log_likelihood_b, num_parameters_b, num_data_points)

print('模型A的BIC:', bic_a.sum())
print('模型B的BIC:', bic_b.sum())

一、利用共轭先验简化计算 #

对于正态似然的场景，对应的共轭先验可以直接给出后验解析解，无需依赖MCMC采样：

• 若标准差σ已知：均值μ的共轭先验为正态分布
• 若标准差σ未知：μ和σ的共轭先验为正态-逆Gamma分布（Normal-Inverse-Gamma）

以未知σ的情况为例，PyMC中可直接定义共轭先验：

with pm.Model() as model_conjugate:
    # 正态-逆Gamma共轭先验：参数根据专家估计调整
    mu = pm.Normal('mu', mu=40, sigma=5)
    sigma = pm.InverseGamma('sigma', alpha=3, beta=20)  # 对应先验均值10
    likelihood = pm.Normal('likelihood', mu=mu, sigma=sigma, observed=data)
    trace_conj = pm.sample(5000, return_inferencedata=True)

二、贝叶斯因子的计算与解读 #

贝叶斯因子是两个模型边缘似然的比值，直接反映数据对模型的支持程度。可通过桥采样在PyMC中计算：

from bridge_sampling import bridge_sampler

# 计算模型A的边缘似然对数
with model_a:
    log_ml_a = bridge_sampler(trace_a).logml

# 计算模型B的边缘似然对数
with model_b:
    log_ml_b = bridge_sampler(trace_b).logml

# 计算贝叶斯因子BF(A/B)
bf_ab = np.exp(log_ml_a - log_ml_b)
print(f"贝叶斯因子BF(A/B): {bf_ab:.2f}")

• 解读规则：BF>3表示数据支持模型A；BF<1/3表示数据支持模型B；1/3<BF<3则数据不足以区分两个模型。

三、LOO值的简化解读 #

LOO（留一法交叉验证）衡量模型的预测性能，az.compare输出的核心指标：

• elpd_loo：模型的期望对数预测密度，值越大预测性能越好
• elpd_diff：模型间elpd_loo的差值，为正说明第一个模型性能更优
• se_diff：差值的标准误，若elpd_diff > 2*se_diff，则差异具有统计显著性

直接打印df_comp_loo即可查看完整对比结果，重点关注上述三列。

四、代码优化建议 #

1. 标准差先验避免用正态分布（可能出现负值），改用HalfNormal或InverseGamma：

sigma = pm.HalfNormal('sigma', sigma=5)  # 仅取正值，更符合标准差的物理意义

2. 采样时添加return_inferencedata=True，让结果兼容ArviZ的所有分析函数：

trace_a = pm.sample(5000, return_inferencedata=True)

3. 修正BIC计算逻辑：手动计算时应使用后验参数均值计算对数似然，而非模型的logp().sum()（后者是先验加似然的对数和，并非边缘似然）：

with model_a:
    # 取后验参数的均值
    mean_post = trace_a.posterior['mean'].mean().item()
    sigma_post = trace_a.posterior['sigma'].mean().item()
    # 计算该参数下的对数似然和
    log_likelihood_a = pm.Normal.dist(mu=mean_post, sigma=sigma_post).logp(data).sum()
bic_a = calculate_bic(log_likelihood_a, 2, len(data))

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