bayesian-regression

贝叶斯回归——Hogg 方法

对于需要可靠的不确定性估计、稳健地处理离群值或异方差误差建模的回归问题，请使用 PyMC 进行贝叶斯回归。本技能采用 Hogg、Bovy 和 Lang（2010）提出的方法：从最简单的模型入手，诊断其失效之处，并升级似然函数以匹配数据实际的噪声结构。

现有的 regression 套件使用 XGBoost 进行非线性表格数据预测。本套件适用于需要具有可信不确定性估计的可解释系数的场景——这类回归常用于论文、监管申报或因果分析。

何时使用本技能

• 需要带有可信区间的系数估计（而不只是点预测）
• 数据存在异方差噪声（方差取决于预测变量）
• 数据存在会导致 OLS/MLE 估计偏差的离群值
• 希望通过有原则的模型比较，比较多个模型设定（Normal 与 Student-t、Poisson）
• 需要通过后验预测检验来验证模型假设
• 希望使用 GLM（逻辑回归、Poisson 回归、Beta 回归）进行完整的贝叶斯推断

何时不使用本技能

• 希望在表格数据上获得尽可能好的预测，而不关心系数 → 使用 regression（XGBoost）套件
• 响应变量是二元的 → 使用 bayesian-ab-testing（转化）或按照下方的 GLM 模式拟合贝叶斯逻辑回归
• 数据具有系数随组别变化的层级/分组结构 → 使用带部分汇聚的 bayesian-regression（参见下方的层级模型部分）

项目结构

<project>/
├── data/                # input parquet/csv
├── src/
│   ├── train.py         # fit models → MLflow log
│   ├── predict.py       # reload idata, posterior predictive
│   └── plots.py         # residuals, PPC, coefficient comparison
├── notebooks/
│   └── demo.py          # marimo walkthrough
└── mlruns/              # MLflow tracking store (gitignored)

Hogg 工作流

步骤 1：OLS 基线

始终从这里开始。OLS 速度快、结果确定，可作为参考。然后检查残差。

import numpy as np

X_design = np.column_stack([np.ones(n), x])
beta_ols = np.linalg.lstsq(X_design, y, rcond=None)[0]
residuals = y - X_design @ beta_ols

残差诊断：

• 绘制残差与各预测变量的关系图 → 扇形表示异方差
• 绘制残差 QQ 图 → 厚尾表示存在离群值
• 如果两者看起来都正常，OLS 就站得住脚。但你仍然没有误差模型的不确定性估计。

步骤 2：贝叶斯 Normal（同方差）

假设与 OLS 相同，但提供完整的后验分布：

import pymc as pm

with pm.Model():
    intercept = pm.Normal("intercept", mu=0, sigma=10)
    slope = pm.Normal("slope", mu=0, sigma=5)
    sigma = pm.HalfNormal("sigma", sigma=5)

    mu = intercept + slope * x
    pm.Normal("y_obs", mu=mu, sigma=sigma, observed=y)

    idata = pm.sample(draws=2000, tune=1000, chains=4)

对预测变量进行中心化。 在拟合前从 x 中减去其均值，可降低截距与斜率之间的后验相关性，使 MCMC 更高效。报告时再还原。

步骤 3：异方差模型

如果残差呈扇形，就将 log(sigma) 建模为 x 的函数：

with pm.Model():
    intercept = pm.Normal("intercept", mu=0, sigma=10)
    slope = pm.Normal("slope", mu=0, sigma=5)

    # Noise model: log(sigma) = gamma0 + gamma1 * x
    gamma0 = pm.Normal("log_sigma_intercept", mu=0, sigma=2)
    gamma1 = pm.Normal("log_sigma_slope", mu=0, sigma=1)
    sigma = pm.math.exp(gamma0 + gamma1 * x)

    mu = intercept + slope * x
    pm.Normal("y_obs", mu=mu, sigma=sigma, observed=y)

这能提供逐点不确定性——数据噪声较大的地方，预测区间更宽；噪声较小的地方，预测区间更窄。OLS 做不到这一点。

步骤 4：稳健模型（Student-t）

Hogg 论文的核心建议：默认使用 Student-t 作为似然函数。 当 nu 较大时，它会退化为 Normal；当 nu 较小时，它则能平稳地降低离群值的权重。

with pm.Model():
    intercept = pm.Normal("intercept", mu=0, sigma=10)
    slope = pm.Normal("slope", mu=0, sigma=5)
    sigma = pm.HalfNormal("sigma", sigma=5)
    nu = pm.Gamma("nu", alpha=2, beta=0.1)  # degrees of freedom

    mu = intercept + slope * x
    pm.StudentT("y_obs", nu=nu, mu=mu, sigma=sigma, observed=y)

解读 nu：

• nu < 5：存在大量离群值污染
• nu ~ 5-30：尾部较厚
• nu > 30：数据基本服从 Gaussian 分布，无需使用 Student-t

步骤 5：通过 LOO-CV 进行模型比较

不要猜哪个模型正确。让预测准确性来决定：

import arviz as az

pm.compute_log_likelihood(idata_homo, model=model_homo)
pm.compute_log_likelihood(idata_hetero, model=model_hetero)
pm.compute_log_likelihood(idata_robust, model=model_robust)

comparison = az.compare({
    "homoscedastic": idata_homo,
    "heteroscedastic": idata_hetero,
    "robust": idata_robust,
}, ic="loo")

elpd_loo 最高的模型具有最佳的样本外预测准确性。weight 给出最优的 stacking 混合权重。

GLM 泛化

只需更改链接函数和似然函数，同一工作流即可扩展到任何 GLM：

响应类型	链接函数	似然函数	PyMC 代码
连续型	identity	Normal / Student-t	pm.Normal("y", mu=mu, sigma=sigma)
二元型	logit	Bernoulli	pm.Bernoulli("y", logit_p=mu)
计数型	log	Poisson	pm.Poisson("y", mu=pm.math.exp(mu))
计数型（过度离散）	log	NegBinomial	pm.NegativeBinomial("y", mu=pm.math.exp(mu), alpha=alpha)
正连续型	log	Gamma	pm.Gamma("y", mu=pm.math.exp(mu), sigma=sigma)
比例 [0,1]	logit	Beta	pm.Beta("y", mu=pm.math.invlogit(mu), kappa=kappa)

诊断工作流相同：拟合、检查 PPC；如果模拟数据与真实数据不匹配，就升级似然函数。

分层 / 部分池化

当数据包含多个组（商店、用户、地区）时，拟合一个分层模型，对系数估计进行部分池化：

with pm.Model(coords={"group": group_names}):
    # Hyperpriors (population-level)
    mu_slope = pm.Normal("mu_slope", mu=0, sigma=5)
    sigma_slope = pm.HalfNormal("sigma_slope", sigma=2)

    # Group-level slopes (partial pooling)
    slope = pm.Normal("slope", mu=mu_slope, sigma=sigma_slope,
                      dims="group")

    # Likelihood
    mu = intercept + slope[group_idx] * x
    pm.Normal("y_obs", mu=mu, sigma=sigma, observed=y)

部分池化会将噪声较大的组估计值向总体均值收缩——数据较多的组保留自己的估计值，数据较少的组则借助总体信息。

诊断——务必检查

1. 收敛性： 所有参数的 R-hat < 1.01，ESS > 400
2. 轨迹图： 链充分混合（呈模糊的毛毛虫状，而不是有趋势或卡住）
3. 后验预测检查： 模拟数据应与真实数据相似
4. 残差图： 模型拟合后不应残留任何模式

MLflow 记录

类型	内容
params	model_type (normal/student_t/hetero)、n_features、draws、tune、chains、seed、centering
metrics	elpd_loo、p_loo、系数（均值 + hdi）、sigma、nu（若为 student-t）、rhat_max、ess_min
tags	data_hash、model_family
artifacts	posterior/idata.nc、plots/{coefficients.png, residuals.png, ppc.png, sigma_x.png, loo_comparison.png}

常见陷阱

1. 默认假设误差服从高斯分布。 将 Student-t 作为默认似然函数，并让数据告诉你高斯分布是否合适（nu >> 30）。
2. 不对预测变量进行中心化。 未中心化的 x 会导致截距与斜率之间出现强后验相关性，从而减慢 MCMC 收敛。务必先进行中心化，再在报告时还原。
3. 忽略异方差。 方差恒定的模型会低估噪声较大处的不确定性，并高估噪声较小处的不确定性。绘制残差与预测变量的关系图。
4. 删除离群值，而不是对其建模。 删除离群值带有主观性，也会丢失信息。Student-t 似然函数会自动降低离群值的权重——让模型来处理它们。
5. 根据系数与 OLS 的相似程度比较模型。 正确的比较方式是使用 LOO-CV 衡量预测准确性，而不是看“哪个模型与 OLS 一致”。整个分析的关键在于，假设不成立时 OLS 并不正确。
6. 不运行后验预测检查。 即使模型的 LOO-CV 表现出色，也仍可能存在 PPC 能揭示的模型设定错误（例如尾部形状不对、遗漏多峰性）。
7. 使用缺乏依据的信息性先验。 从弱信息先验开始（系数使用 Normal(0, 10)，尺度使用 HalfNormal(5)）。只有在掌握真实领域知识时才收紧先验。

示例

参见 demo.py（marimo notebook）。它会生成包含预设异方差和离群值的合成数据，拟合 OLS 和三个贝叶斯模型（Normal、异方差、Student-t），通过 LOO-CV 进行比较，并展示后验预测检查。运行命令：

marimo edit --sandbox demo.py

参考资料

• Hogg、Bovy 和 Lang（2010），《数据分析秘诀：拟合数据模型》（arXiv:1008.4686）
• McElreath（2020），统计学思维重构，第 4–5 章（回归）和第 13–14 章（分层模型）
• Gelman 等人（2013），贝叶斯数据分析，第 14 章（GLM）