PyMC 贝叶斯分位数回归：3 步让安全库存不再只靠均值

PyMC 贝叶斯分位数回归：3 步让安全库存不再只靠均值

【免费下载链接】pymcBayesian Modeling and Probabilistic Programming in Python项目地址: https://gitcode.com/GitHub_Trending/py/pymc

备多少货才不缺货？均值预测永远答不准：它只给出“平均卖多少”，而你要的是“最坏卖多少”。PyMC 的贝叶斯分位数回归把 90%、95% 分位数直接当成回归目标，用 MCMC 后验把“卖爆”风险算成一条可解释的曲线。

分位数回归到底在解决什么问题

• 金融风控：真正危险的是尾部损失，95% 分位数才是预警线，平均损失只是参考。
• 供应链：按平均需求备货，大促日必然缺货，90% 分位数才配得上“安全库存”四个字。
• 用户留存：平均使用时长分不清大 R 和路人，分位数能把用户切进不同运营策略的档位。

直觉上，分位数就是“按位置取数”。班级平均分 70，但按排名看，75 分位可能是 82 分——均值压平了分布，分位数保留了结构。分位数回归做的事，就是把“每个 x 对应一条平均线”，换成“每个 x 对应一组指定位置的分位线”。

PyMC 里的贝叶斯做法：先验、似然、后验

三步走，每步一句话：

• 先验：把业务经验写成参数范围，比如斜率不会超过 ±5。
• 似然：描述数据如何从参数中生成，这里选不对称拉普拉斯分布。
• 后验：MCMC 采样出所有参数的联合分布，不确定性量化是副产品。

PyMC 把分布、采样器和 ArviZ 诊断组织在同一条流水线里：

似然选不对称拉普拉斯分布（AsymmetricLaplace），它把误差按分位数加权：

$$f(y\mid\mu,\sigma,q)\propto\exp\left(-\frac{|y-\mu|}{\sigma}\left[q,\mathbb{1}(y\ge\mu)+(1-q),\mathbb{1}(y<\mu)\right]\right)$$

q=0.9 时，“预测偏低”的代价是“预测偏高”的 9 倍，μ 自然被推到 90% 分位。建模写法很小：

import numpy as np import pandas as pd import pymc as pm import arviz as az x_data = np.linspace(0, 10, 100) y = 1 + 1.2 * x_data + np.random.normal(0, 1, 100) with pm.Model() as model: beta = pm.Normal("beta", mu=0, sigma=5) # 宽先验：先让数据说话 sigma = pm.HalfNormal("sigma", sigma=2) # q 决定位置参数最终收敛到哪个分位数 pm.AsymmetricLaplace("y_obs", mu=1.0 + beta * x_data, b=sigma, q=0.9, observed=y)

从数据到诊断：PyMC 分位数回归实操流程

先造一批异方差模拟数据，方差随 x 变大，正是均值回归吃亏的场景：

np.random.seed(42) n = 300 x = np.linspace(0, 10, n) y = 2 + 1.2 * x + np.random.normal(0, 0.8 * (x / 10 + 0.5), n)

同一份数据上建 90% 分位数模型并采样：

with pm.Model() as model: beta = pm.Normal("beta", mu=0, sigma=5) sigma = pm.HalfNormal("sigma", sigma=2) mu = 1.0 + beta * x pm.AsymmetricLaplace("y_obs", mu=mu, b=sigma, q=0.9, observed=y) idata = pm.sample(2000, cores=2, tune=1000)

诊断只盯两件事：R-hat 是否收敛、后验预测区间覆盖率是否接近名义分位：

print(idata.sample_stats.r_hat) # 全部 < 1.01 才算收敛 az.plot_forest(idata) # 森林图：可信区间 + r_hat 一览 ppc = pm.sample_posterior_predictive(idata) lo, hi = np.percentile(ppc.posterior_predictive.y_obs, [5, 95], axis=(0, 1)) print(float(np.mean((y > lo) & (y < hi)))) # ≈0.9 说明模型可信

森林图里，94% 可信区间窄、r_hat 贴近 1，说明截距和斜率都估稳了：

多分位数联合建模与贝叶斯不确定性量化

分开跑三次，得到的是三组互不相干的参数：无法互相比较，也无法保证 10% 线始终在 90% 线下方。联合建模让三条分位线共享同一批数据，后验里天然带着分位间的相关结构。

qs = [0.1, 0.5, 0.9] with pm.Model() as model: beta = pm.Normal("beta", mu=0, sigma=5, shape=len(qs)) sigma = pm.HalfNormal("sigma", sigma=2, shape=len(qs)) for i, q in enumerate(qs): pm.AsymmetricLaplace(f"y_q{int(q*100)}", mu=1.0 + beta[i] * x, b=sigma[i], q=q, observed=y) idata = pm.sample(2000, cores=2) post = idata.posterior["beta"].mean(("chain", "draw"))

把post[0]和post[2]画出来就是 10% 与 90% 分位数曲线。图上看点：50% 线贴近均值回归线；10% 与 90% 线的间距随 x 拉大，说明条件分布的离散度在增长；间距本身就是“此处预测更不可靠”的可视化。

💡 两个可以照搬的落地模板

场景 A：客户 LTV 的 90% 分位数预测，重点在宽先验与显式提取分位数：

df = pd.read_csv("customer_ltv.csv") # 字段：recency, frequency, ltv X = df[["recency", "frequency"]].values with pm.Model() as ltv_model: alpha = pm.Normal("alpha", mu=0, sigma=10) beta = pm.Normal("beta", mu=0, sigma=10, shape=2) sigma = pm.HalfNormal("sigma", sigma=10) pm.AsymmetricLaplace("ltv", mu=alpha + X @ beta, b=sigma, q=0.9, observed=df["ltv"].values) idata_ltv = pm.sample(1500) p90 = idata_ltv.posterior["ltv"] # 逐客户的 90% LTV 后验，直接分群

上线前最需要检查：拿历史数据回测，90% 预测值高于真实 LTV 的客户比例是否约等于 90%。

场景 B：零售需求的 95% 上界用于补货点，分层结构给每个门店一个截距：

df = pd.read_csv("demand.csv") # 字段：store, day, demand store_idx = pd.factorize(df["store"])[0] with pm.Model() as inv_model: store_mu = pm.Normal("store_mu", mu=50, sigma=20, shape=df["store"].nunique()) beta = pm.Normal("beta", mu=0, sigma=5) sigma = pm.HalfNormal("sigma", sigma=10) pm.AsymmetricLaplace("demand", mu=store_mu[store_idx] + beta * df["day"], b=sigma, q=0.95, observed=df["demand"].values) idata_inv = pm.sample(1500) p95 = idata_inv.posterior["demand"].mean(("chain", "draw")) # 直接进补货点公式

上线前最需要检查：95% 区间对历史真实需求的覆盖率，以及门店截距有没有被先验过度收缩。

⚠️ 常见坑与下一步

• 先验过宽：Normal(sigma=100) 会让后验又平又钝，MCMC 走得慢、R-hat 容易破 1.01；业务上能约束的，一律写紧先验。
• 分位数交叉：各分位数斜率独立时，大 x 处 10% 线可能越过 90% 线；用共享结构或对斜率加单调约束。
• 异方差未处理：方差随 x 变化时，单一 sigma 必然顾此失彼；先让 sigma 成为协变量的函数，再谈拟合好坏。
• q 太贴边：q=0.999 的似然极度偏斜，自相关飙升；先用 0.99 验证，再谈极端分位。

进阶方向：非线性分位数回归，把 beta·x 换成样条基展开，PyMC 里只是多一个点积；时序分位数，给分位函数加高斯过程先验，捕捉分位数随时间的漂移。

API 速查（AsymmetricLaplace）：

延伸阅读：PyMC 官方文档、概率分布指南、连续分布单元测试。

【免费下载链接】pymcBayesian Modeling and Probabilistic Programming in Python项目地址: https://gitcode.com/GitHub_Trending/py/pymc