PyMC5贝叶斯回归模型：如何获取AIC|BIC|LOO用于模型比较

一、先算LOO（留一交叉验证）评分 #

LOO是贝叶斯模型对比里更靠谱的指标，PyMC5有现成工具，步骤简单：

1. 采样后直接计算对数似然和LOO：

with basic_model:
    # 已经采样过的话跳过这行
    idata = pm.sample(1000)
    # 计算对数似然
    idata = pm.compute_log_likelihood(idata)
    # 算出LOO
    loo_result = pm.loo(idata)
    print(loo_result)

输出里的LOO数值越小（越负），说明模型拟合效果越好，用来对比不同模型很方便。

二、计算AIC和BIC #

你之前的报错是因为处理标量参数时用了len()，标量没有长度属性，换成np.size()就能解决，完整步骤：

1. 获取MAP（最大后验概率）估计值 #

with basic_model:
    map_est = pm.find_MAP()

2. 统计总参数数量k #

遍历MAP结果里的每个参数，用np.size()统计每个参数的元素个数（不管是标量还是数组都能处理）：

import numpy as np

k = sum(np.size(v) for v in map_est.values())

3. 计算对数似然值 #

编译模型的logp函数，传入MAP估计值得到对数似然：

f_logp = basic_model.compile_logp()
logp_value = f_logp(map_est)

4. 套公式计算AIC和BIC #

n = len(y_array)  # 你的观测样本量
aic = 2 * k - 2 * logp_value
bic = k * np.log(n) - 2 * logp_value

print(f"AIC: {aic:.2f}")
print(f"BIC: {bic:.2f}")

三、整合后的完整代码 #

把这些步骤加到你的代码里：

import pymc as pm
import numpy as np

# 假设你的a_array、b_array、c_array、y_array已经定义好
# ... 这里放你的数据定义代码 ...

# 定义模型
basic_model = pm.Model()

with basic_model:
    cost = pm.Normal("cost", mu=0, sigma=2)
    a = pm.Normal('a_array', mu=0, sigma=2)
    b = pm.Normal('b_array', mu=0, sigma=2)
    c = pm.HalfNormal('c_array', sigma=2)
    sigma = pm.HalfNormal("sigma", sigma=1)

    mu = cost + a * a_array + b * b_array + c * c_array
    Y_obs = pm.Normal("Y_obs", mu=mu, sigma=sigma, observed=y_array)

# 计算LOO
with basic_model:
    idata = pm.sample(1000)
    idata = pm.compute_log_likelihood(idata)
    loo_result = pm.loo(idata)
    print("LOO结果：")
    print(loo_result)

# 计算AIC和BIC
with basic_model:
    map_est = pm.find_MAP()
    k = sum(np.size(v) for v in map_est.values())
    f_logp = basic_model.compile_logp()
    logp_value = f_logp(map_est)
    n = len(y_array)
    aic = 2 * k - 2 * logp_value
    bic = k * np.log(n) - 2 * logp_value
    print(f"\nAIC: {aic:.2f}")
    print(f"BIC: {bic:.2f}")

要点提醒 #

• 优先用LOO做模型对比，它考虑了后验的不确定性，比AIC/BIC更适合贝叶斯框架。
• pm.find_MAP()默认用BFGS优化器，你的回归模型用默认设置就行，复杂模型可能需要调整优化参数。
• 之前的len() of unsized object错误，就是因为cost这种标量参数用len()获取长度，换成np.size()就解决了。

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