均值骗人：PyMC 贝叶斯分位数回归抓需求上界

均值骗人：PyMC 贝叶斯分位数回归抓需求上界

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

大促前主管拍板：备货按九成九不缺货来定。你算了日均 850 件就照它下单，结果第三天断货——平均数压不住促销期抬起来的尾部。要的是一条九成概率不超它的线，也就是条件分位数。用 PyMC 贝叶斯分位数回归，就是直接把这条上界建模出来。读完你能写出 AsymmetricLaplace 模型、看懂收敛诊断，再把后验区间翻译成能下单的备货量。

一个帮你建立直觉的比喻

先看个熟悉的东西：体检报告上的某项指标，医生不会只盯着"人群平均"，而是看你的值落在第几分位——低于 P10 要警惕，高于 P90 也要警惕。均值只告诉你人群中心在哪，分位数告诉你整条分布长什么样、尾巴有多长。

天气预报更直白。它说"明天降雨 5 毫米"对你没用，你要的是"会不会大到淋湿"这种阈值判断，而阈值判断靠的就是分位。

回归也一样。传统线性回归给的是条件均值，也就是 Y 的中心点。分位数回归给的是整条条件分布的轮廓：你可以同时描出 P10、P50、P90 三条线，把典型值和极端值一起框住。说白了，非正态分布回归想做的事就是不逼数据服从正态，而是承认分布有形状，然后直接去描这个形状。你可以把它理解为：均值是拍一张正面照，分位数是绕着物体转一圈。

30 秒跑通第一个模型

下面这段能直接跑：造一条带异方差噪声的线，再用AsymmetricLaplace把 90% 分位当似然拟合出来。数据怎么造的一行带过，重点看分布定义那几行。

import numpy as np import pymc as pm rng = np.random.default_rng(0) x = rng.uniform(0, 10, 300) # 自变量，比如促销力度 y = 2.0 + 0.8 * x + rng.normal(0, 0.5 * (x / 10 + 0.3), 300) # 模拟数据：噪声随 x 变大 with pm.Model() as m: beta0 = pm.Normal("beta0", 0, 10) # 截距先验 beta1 = pm.Normal("beta1", 0, 10) # 斜率先验 sigma = pm.HalfNormal("sigma", 2) # 尺度参数先验 mu = beta0 + beta1 * x # 线性预测：条件分位数函数 pm.AsymmetricLaplace("y_obs", mu=mu, b=sigma, q=0.9, observed=y) # q=0.9 建模 90% 分位 idata = pm.sample(1000, tune=500, progressbar=False) # MCMC 采样

三个参数各管一件事：mu就是你要预测的那条分位线本身，业务上是"九成不超它"的需求量；sigma是分布的离散程度，越大说明同一个 x 下结果越散；q决定你盯哪条线，q=0.9就是 90% 分位，想看中位数就填 0.5。这里用 AsymmetricLaplace 当似然，是因为它天然把 q 编码进损失里，q 越大越重罚"预测偏低"，正好匹配"宁可多备也别缺货"的心态。这个分布的实现就放在 continuous.py 里。

PyMC 会把上面这些节点自动连成一张概率图模型，随机变量是自由节点，派生量挂在它们下游。

怎么判断采样结果靠谱

采样完别急着用，先看三个信号。

第一是 R-hat，衡量各条链有没有混到同一处，落在 1.01 以内算稳，超过 1.1 就该重跑。第二是 ESS（有效样本量），它不是越大越好的装饰，而是"你这几千个样本里有多少个真正独立"，太低说明样本高度自相关，后验区间不可信。第三是后验预测图：让模型重新生成一批假数据，叠到真实数据上看形状和尾部对不对得上，对得上才说明模型真的描述了数据。

一行命令把前两样调出来：

import arviz as az az.plot_trace(idata, var_names=["beta0", "beta1", "sigma"])

上图左半是每个参数的可信区间（粗横线是后验主体），右半是 R-hat 点——都贴着 1，说明这条链收敛了。读图比调参重要：先看右半收没收敛，再看左半区间宽不宽。

从一条线到三条线：多分位数建模

业务里你很少只关心一条线。客服排班看的是中位数，风险盯的是 P95，清仓促销要看 P10。所以条件分位数建模常常一次画三条，好把"典型、偏高、偏低"一起摆出来。

技巧就一招：用shape把参数拉成一维向量，每条分位线配一套自己的beta和sigma，再循环绑定似然。

qs = [0.1, 0.5, 0.9] # 同时看 10%、50%、90% 三条线 with pm.Model() as mq: beta0 = pm.Normal("beta0", 0, 10, shape=len(qs)) # 每条线一个截距 beta1 = pm.Normal("beta1", 0, 10, shape=len(qs)) # 每条线一个斜率 sigma = pm.HalfNormal("sigma", 2, shape=len(qs)) mu = beta0 + beta1[None, :] * x[:, None] # 广播成 (n_obs, n_quantiles) for i, q in enumerate(qs): # 循环给每条线配似然 pm.AsymmetricLaplace(f"y_{q}", mu=mu[:, i], b=sigma[i], q=q, observed=y)

这里有个反直觉的点：三条线不是各画各的，它们共享同一套数据，MCMC 会把它们之间的相关性一并学出来。画出来你会看到三条线像扇面一样张开——x 小时彼此贴着，x 大时上下分位离中位越来越远，把"越往右越不确定"这件事直接描在了图上。

一个落地场景：电商补货的需求上界

回到开头那个断货的坑。补货决策的核心诉求是需求上限，不是均值。建模思路分三步。

特征怎么选：促销开关、折扣力度、星期几、天气，这些直接影响当天需求；但别急着塞"上周销量"进去，它和当前需求高度相关，容易把共线性搅进后验。

为什么取q=0.95而不是0.90：因为缺货和滞销的成本不对称。大促断货会丢流量、招差评，多压一点库存只是占用资金，代价低得多，所以上界要往更保守的那端取。

怎么把后验翻译成备货量：别拿点估计下单，取后验的一个分位当安全备货量。

X = feats[["promo", "discount", "dow"]] # 特征：促销 / 折扣 / 星期 with pm.Model() as om: beta = pm.Normal("beta", 0, 5, shape=X.shape[1]) # 每个特征一个系数 alpha = pm.Normal("alpha", 0, 10) sigma = pm.HalfNormal("sigma", 100) mu = alpha + X.dot(beta) # 线性预测：95% 需求上界 pm.AsymmetricLaplace("demand", mu=mu, b=sigma, q=0.95, observed=demand) idata = pm.sample()

举个具体的数：模型给出某天需求上界的后验，95% 可信区间是 [800, 1200]，中位数 900。如果按均值 900 备货，等于赌"需求不会超过 900"，而它超过 900 的概率不小，粗算有 25% 左右会断货；反过来直接按 1200 备，又可能白压三成库存成本。正确做法是取后验的 95% 分位（比如 1150）当安全备货量，把断货风险显式压到你愿意接受的 5%。这就是概率化预测比"甩一个数"值钱的地方。

三个常见坑

第一，q 贴边。现象是q设成 0.99 甚至 1.0 时参数检查直接报错，或后验剧烈漂移。原因是 q 必须严格落在 0 到 1 之间，PyMC 内部把它换算成对称参数 kappa，q 越接近 0 或 1，kappa 越极端，分布退成单边、梯度爆炸。避免办法是业务分位留点余地，0.95、0.99 都行，但别用 1.0，真要够尾部就加大 tune。

第二，把 q 当置信区间。现象是有人直接拿q=0.95的 mu 后验当"95% 置信区间"去汇报。原因是 q 描述的是数据分布的分位，即需求落在它以下的概率，而置信区间讲的是参数估计本身的不确定性，两码事。避免办法是分位线只用来报"预测需求的概率上界"，参数的不确定性另看 trace 图上的可信带，别混着讲。

第三，shape 对齐出错。现象是多分位数一跑就报形状不匹配，mu 广播对不上。原因是beta0、beta1用了shape拉成向量，但 mu 里索引的维度没对齐，或循环里取错了轴。避免办法是先在小数据集上把 shape 打印出来核对，确认 mu 是 (n_obs, n_quantiles) 再喂给似然。

接下来可以做什么

三条能马上做的延伸。第一，嵌进 A/B 测试：把分位数回归当实验分析的一步，直接比较新旧策略在不同分位上的差异，而不只看均值提升。第二，试非线性：x 和分位线关系弯曲时，把线性项换成样条，或用分层模型同时给多个门店、区域各画一套分位线。第三，接时间序列：需求有周期和趋势时，把 AsymmetricLaplace 似然挂到状态空间模型上，让分位线随时间走。

想继续深挖，两个地方值得一读：概率分布指南 里对不对称拉普拉斯的展开，和 GLM 线性回归教程 用来对比传统回归和分位数回归的差别。

PyMC 把分布、采样、诊断拆成了独立模块，你以后想换分布、换采样器，都是在这张图里挪一下的事。把这套流程跑通，你就已经超过了大多数还只盯着 OLS 残差的人。

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