When does Bayes win? — GBM vs a hierarchical pooling model¶
Machine Learning in Operations Research · the cold-start crossover¶
The first two notebooks established that quantile gradient boosting is an excellent, fast distributional forecaster for this inventory problem — it won on the best-sellers and scaled to the full catalogue. But GBM leans on data: it learns each product's behaviour from that product's history. When history is thin — a newly listed item, a rare slow mover — there is little to learn from.
That is exactly where a Bayesian hierarchical model with product random effects should shine: with little data, partial pooling shrinks each product toward the population, and priors substitute for missing observations. This notebook races the two head-to-head across data-richness strata to see whether, and where, the Bayesian model overtakes GBM.
The two contenders (same leak-safe features, same 28-day holdout, scored by realized newsvendor £):
- Quantile GBM — one LightGBM per quantile (notebooks 1–2), product identity via a native category.
- Bayesian hierarchical NegBin GLMM — $\text{demand}\sim\text{NegBin}(\mu_{pt},\phi)$, $\log\mu_{pt}=a_p+x_{pt}^\top\beta$, with product random intercepts $a_p\sim\mathcal N(\mu_a,\sigma_a)$. The $a_p$ are where partial pooling happens. Fit with NUTS in NumPyro.
In [1]:
import os, time
os.environ["KMP_DUPLICATE_LIB_OK"] = "TRUE"
import numpy as np, pandas as pd, matplotlib.pyplot as plt
from pathlib import Path
BLUE, RED, GREEN, ORANGE, GREY, PURP = "#2b6cb0","#c53030","#2f855a","#dd6b20","#a0aec0","#6b46c1"
here = Path.cwd(); DATA = next(p for p in [here/"data", here.parent/"data"] if (p/"features_full.parquet").exists())
feat = pd.read_parquet(DATA/"features_full.parquet"); feat["date"] = pd.to_datetime(feat["date"])
MARGIN_RATE, HOLDING = 0.5, 0.10
QLEV = [0.1, 0.25, 0.5, 0.7, 0.8, 0.9, 0.95, 0.99]
price = feat.groupby("StockCode", observed=True)["price"].first(); cu_prod = MARGIN_RATE*price; fractile = cu_prod/(cu_prod+HOLDING)
# --- stratified sample of ~150 products spanning history-richness ---
cut = feat["date"].max() - pd.Timedelta(days=28)
train_all = feat[feat["date"] <= cut]
active = train_all[train_all["demand"] > 0].groupby("StockCode", observed=True).size() # nonzero days in train
active = active[active >= 5]
rng = np.random.default_rng(7)
terc = pd.qcut(active, 3, labels=["scarce", "medium", "rich"])
picks = []
for lab in ["scarce", "medium", "rich"]:
pool = active[terc == lab].index.to_numpy()
picks += list(rng.choice(pool, size=min(50, len(pool)), replace=False))
picks = pd.Index(picks)
strat = terc.reindex(picks)
print(f"{len(picks)} products | active-days by stratum: "
f"scarce {active.reindex(strat[strat=='scarce'].index).median():.0f}, "
f"medium {active.reindex(strat[strat=='medium'].index).median():.0f}, "
f"rich {active.reindex(strat[strat=='rich'].index).median():.0f} (median nonzero days)")
150 products | active-days by stratum: scarce 20, medium 62, rich 183 (median nonzero days)
The data — a subset spanning the richness spectrum¶
We keep the same leak-safe features as before, restricted to these products, and drop warm-up rows.
In [2]:
USE = ["is_saturday", "sin_year", "cos_year", "days_to_xmas", "roll_mean_28", "roll_zero_28", "lag_7"]
d = feat[feat.StockCode.isin(picks)].copy()
d["logrm"] = np.log1p(d["roll_mean_28"]); d["loglag7"] = np.log1p(d["lag_7"])
USEc = ["is_saturday","sin_year","cos_year","days_to_xmas","logrm","roll_zero_28","loglag7"]
d = d.dropna(subset=USEc + ["demand"]).copy()
tr, te = d[d["date"] <= cut], d[d["date"] > cut]
te = te.assign(stratum=te["StockCode"].map(strat).astype(str))
print(f"train rows {len(tr):,} | holdout rows {len(te):,} | holdout products {te.StockCode.nunique()}")
def order_at(P, fr, qs):
q = np.array(sorted(qs)); V = P[sorted(qs)].to_numpy()
return np.array([np.interp(fr[i], q, V[i]) for i in range(len(fr))])
def pinball(y, q, tau): e = y-q; return float(np.mean(np.maximum(tau*e, (tau-1)*e)))
def costs_by_stratum(P, te):
fr = fractile.reindex(te["StockCode"]).to_numpy()
order = order_at(P.loc[te.index], fr, QLEV)
D = te["demand"].to_numpy(float); CU = cu_prod.reindex(te["StockCode"]).to_numpy()
row_cost = CU*np.maximum(D-order, 0) + HOLDING*np.maximum(order-D, 0)
out = pd.DataFrame({"stratum": te["stratum"].to_numpy(), "cost": row_cost, "demand": D})
g = out.groupby("stratum").apply(lambda x: pd.Series({"cost": x["cost"].sum(),
"cost_per_demand": x["cost"].sum()/max(x["demand"].sum(), 1)}), include_groups=False)
return g
train rows 104,550 | holdout rows 4,200 | holdout products 150
Fit the two models¶
In [3]:
import lightgbm as lgb
t = time.time()
cat = ["StockCode", "sb_class"]
Xtr = tr[USEc + cat].copy(); Xte = te[USEc + cat].copy()
for c in cat: Xtr[c] = Xtr[c].astype("category"); Xte[c] = Xte[c].astype("category")
cv = tr["date"].max() - pd.Timedelta(days=28); fitm = tr[tr.date <= cv]; valm = tr[tr.date > cv]
Xf, Xv = fitm[USEc+cat].copy(), valm[USEc+cat].copy()
for c in cat: Xf[c]=Xf[c].astype("category"); Xv[c]=Xv[c].astype("category")
gq = {}
for tau in QLEV:
m = lgb.LGBMRegressor(objective="quantile", alpha=tau, n_estimators=300, learning_rate=0.05,
num_leaves=31, min_child_samples=20, verbosity=-1)
m.fit(Xf, fitm["demand"], eval_set=[(Xv, valm["demand"])], eval_metric="quantile",
callbacks=[lgb.early_stopping(40, verbose=False), lgb.log_evaluation(0)]); gq[tau] = m.predict(Xte)
Pg = np.clip(np.column_stack([gq[t_] for t_ in QLEV]), 0, None); Pg.sort(axis=1)
Pg = pd.DataFrame(Pg, columns=QLEV, index=te.index)
print(f"GBM fit in {time.time()-t:.1f}s")
C:\Users\user\anaconda3\envs\pymc-env\Lib\site-packages\lightgbm\sklearn.py:1106: LGBMDeprecationWarning: The argument 'eval_set' is deprecated, use 'eval_X' and 'eval_y' instead. eval_set = _validate_eval_set_Xy(eval_set=eval_set, eval_X=eval_X, eval_y=eval_y) C:\Users\user\anaconda3\envs\pymc-env\Lib\site-packages\lightgbm\sklearn.py:1106: LGBMDeprecationWarning: The argument 'eval_set' is deprecated, use 'eval_X' and 'eval_y' instead. eval_set = _validate_eval_set_Xy(eval_set=eval_set, eval_X=eval_X, eval_y=eval_y)
C:\Users\user\anaconda3\envs\pymc-env\Lib\site-packages\lightgbm\sklearn.py:1106: LGBMDeprecationWarning: The argument 'eval_set' is deprecated, use 'eval_X' and 'eval_y' instead. eval_set = _validate_eval_set_Xy(eval_set=eval_set, eval_X=eval_X, eval_y=eval_y)
C:\Users\user\anaconda3\envs\pymc-env\Lib\site-packages\lightgbm\sklearn.py:1106: LGBMDeprecationWarning: The argument 'eval_set' is deprecated, use 'eval_X' and 'eval_y' instead. eval_set = _validate_eval_set_Xy(eval_set=eval_set, eval_X=eval_X, eval_y=eval_y)
C:\Users\user\anaconda3\envs\pymc-env\Lib\site-packages\lightgbm\sklearn.py:1106: LGBMDeprecationWarning: The argument 'eval_set' is deprecated, use 'eval_X' and 'eval_y' instead. eval_set = _validate_eval_set_Xy(eval_set=eval_set, eval_X=eval_X, eval_y=eval_y)
C:\Users\user\anaconda3\envs\pymc-env\Lib\site-packages\lightgbm\sklearn.py:1106: LGBMDeprecationWarning: The argument 'eval_set' is deprecated, use 'eval_X' and 'eval_y' instead. eval_set = _validate_eval_set_Xy(eval_set=eval_set, eval_X=eval_X, eval_y=eval_y)
C:\Users\user\anaconda3\envs\pymc-env\Lib\site-packages\lightgbm\sklearn.py:1106: LGBMDeprecationWarning: The argument 'eval_set' is deprecated, use 'eval_X' and 'eval_y' instead. eval_set = _validate_eval_set_Xy(eval_set=eval_set, eval_X=eval_X, eval_y=eval_y)
C:\Users\user\anaconda3\envs\pymc-env\Lib\site-packages\lightgbm\sklearn.py:1106: LGBMDeprecationWarning: The argument 'eval_set' is deprecated, use 'eval_X' and 'eval_y' instead. eval_set = _validate_eval_set_Xy(eval_set=eval_set, eval_X=eval_X, eval_y=eval_y)
GBM fit in 5.3s
In [4]:
import jax, jax.numpy as jnp, numpyro, numpyro.distributions as dist
from numpyro.infer import MCMC, NUTS, Predictive
pid = {p: i for i, p in enumerate(picks)}
mu_, sd_ = tr[USEc].to_numpy(float).mean(0), tr[USEc].to_numpy(float).std(0)+1e-6
Ztr = (tr[USEc].to_numpy(float)-mu_)/sd_; Zte = (te[USEc].to_numpy(float)-mu_)/sd_
ptr = tr["StockCode"].map(pid).to_numpy().astype(np.int32); pte = te["StockCode"].map(pid).to_numpy().astype(np.int32)
def model(prod, X, y=None):
mu_a = numpyro.sample("mu_a", dist.Normal(0., 3.)); sig_a = numpyro.sample("sig_a", dist.HalfNormal(3.))
with numpyro.plate("p", len(picks)): a = numpyro.sample("a", dist.Normal(mu_a, sig_a))
beta = numpyro.sample("beta", dist.Normal(0., 1.).expand([X.shape[1]])); conc = numpyro.sample("conc", dist.HalfNormal(5.))
numpyro.sample("obs", dist.NegativeBinomial2(jnp.exp(a[prod] + X @ beta), conc), obs=y)
t = time.time()
mcmc = MCMC(NUTS(model), num_warmup=350, num_samples=350, num_chains=1, progress_bar=False)
mcmc.run(jax.random.PRNGKey(0), jnp.asarray(ptr, jnp.int32), jnp.asarray(Ztr), y=jnp.asarray(tr["demand"].to_numpy(float)))
pred = np.asarray(Predictive(model, mcmc.get_samples())(jax.random.PRNGKey(1), jnp.asarray(pte, jnp.int32), jnp.asarray(Zte))["obs"])
Pb = pd.DataFrame({q: np.quantile(pred, q, axis=0) for q in QLEV}, index=te.index)
print(f"Bayesian NUTS fit in {time.time()-t:.0f}s | {len(picks)} product random effects")
Bayesian NUTS fit in 111s | 150 product random effects
The crossover — realized cost by data-richness¶
We score both models' realized newsvendor cost within each stratum, normalized by demand so strata are comparable. The hypothesis: GBM leads where data is rich; the Bayesian model closes the gap — or overtakes — where data is scarce.
In [5]:
g_gbm = costs_by_stratum(Pg, te); g_bay = costs_by_stratum(Pb, te)
comp = pd.DataFrame({"GBM": g_gbm["cost_per_demand"], "Bayes-RE": g_bay["cost_per_demand"]}).reindex(["scarce","medium","rich"])
comp["Bayes/GBM"] = (comp["Bayes-RE"]/comp["GBM"]).round(3)
print("realized cost per unit demand, by stratum:"); print(comp.round(3).to_string())
print("\noverall pinball q0.9: GBM %.2f | Bayes %.2f" % (pinball(te.demand.values, Pg[0.9].values, .9), pinball(te.demand.values, Pb[0.9].values, .9)))
fig, ax = plt.subplots(1, 2, figsize=(13, 4.5))
x = np.arange(3); w = 0.38
ax[0].bar(x-w/2, comp["GBM"], w, color=BLUE, label="GBM")
ax[0].bar(x+w/2, comp["Bayes-RE"], w, color=GREEN, label="Bayes + product RE")
ax[0].set_xticks(x); ax[0].set_xticklabels(comp.index); ax[0].set_ylabel("realized cost per unit demand (£)")
ax[0].set_title("Cost by data-richness (lower = better)"); ax[0].legend(fontsize=8)
ax[1].plot(comp.index, comp["Bayes/GBM"], "o-", color=PURP, lw=2); ax[1].axhline(1.0, color=GREY, ls="--", label="parity")
ax[1].set_ylabel("Bayes cost ÷ GBM cost"); ax[1].set_title("Bayesian advantage grows as data thins"); ax[1].legend(fontsize=8)
for i, v in enumerate(comp["Bayes/GBM"]): ax[1].annotate(f"{v:.2f}", (i, v), xytext=(0,6), textcoords="offset points", ha="center", fontsize=9)
fig.tight_layout(); plt.show()
realized cost per unit demand, by stratum:
GBM Bayes-RE Bayes/GBM
stratum
scarce 1.010 0.713 0.706
medium 0.340 0.635 1.870
rich 0.301 0.664 2.205
overall pinball q0.9: GBM 3.14 | Bayes 4.70
Verdict¶
Read the right-hand panel: the ratio Bayes ÷ GBM falls as products get scarcer — the Bayesian model's relative position improves exactly where history is thin, because partial pooling and priors do the work that data cannot. On the data-rich products, GBM's flexibility wins; on the scarce tail, the hierarchical model is competitive or better.
The practical takeaway — a hybrid. Neither model dominates everywhere, so the right production design routes by data-richness:
- data-rich products → quantile GBM (fast, flexible, best where signal is plentiful);
- sparse / cold-start products → the Bayesian hierarchical model (pooling + priors substitute for missing history).
That plays to both halves of the toolkit — the ML forecaster of notebooks 1–2 and the Bayesian hierarchical modelling from the broader portfolio — and it is the honest answer to "which model is best?": it depends on how much data the product has.
(A DeepAR-style probabilistic RNN was also benchmarked separately and landed between the two on the data-rich set, while requiring per-series scaling to train at all — competitive, but not a clear win here, consistent with GBMs' edge on tabular retail.)