Predict-then-Optimize at scale — the full catalogue¶
Machine Learning in Operations Research · the long tail, 4,862 products¶
This is the companion to the top-200 notebook (predict_then_optimize_uci.ipynb). Same five-act
pipeline — clean → distributional forecast → newsvendor decision → backtest — but run on the entire
UK catalogue rather than the best-sellers. The point is what changes when you leave the easy,
high-volume products and take on the long tail of slow movers, which is where real inventory
management lives and where distributional forecasting matters most.
Where the method is identical, we keep it brief and point back to the top-200 notebook for the full derivations (newsvendor rule, pinball loss, leakage discipline, LP + shadow prices). The focus here is the difference at scale.
In [1]:
import os, re, 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/"online_retail_II.parquet").exists())
raw = pd.read_parquet(DATA/"online_retail_II.parquet")
print(f"{len(raw):,} raw transactions")
1,067,371 raw transactions
1 · Clean → the full demand panel¶
Identical cleaning to the top-200 notebook (UK · drop cancellations & non-positive quantity · real product codes · positive price), but we keep every product, not the top 200.
In [2]:
def clean(df):
df = df[(df.Country=="United Kingdom") & (~df.Invoice.str.startswith("C")) & (df.Quantity>0) & (df.Price>0)].copy()
df = df[df.StockCode.str.match(r"^\d{5}[A-Za-z]?$")]
df["date"] = df.InvoiceDate.dt.normalize()
return df
c = clean(raw)
daily = c.groupby(["StockCode","date"], as_index=False)["Quantity"].sum()
prods = daily.StockCode.unique(); all_days = pd.date_range(c.date.min(), c.date.max(), freq="D")
panel = (daily.set_index(["StockCode","date"])["Quantity"]
.reindex(pd.MultiIndex.from_product([prods, all_days], names=["StockCode","date"]), fill_value=0)
.rename("demand").reset_index())
panel["description"] = panel.StockCode.map(c.groupby("StockCode")["Description"].agg(lambda s: s.mode().iloc[0]))
panel["price"] = panel.StockCode.map(c.groupby("StockCode")["Price"].median())
print(f"panel: {len(panel):,} rows | {panel.StockCode.nunique():,} products x {len(all_days)} days | "
f"zero-share {(panel.demand==0).mean():.0%} | median price GBP {panel.price.median():.2f}")
panel: 3,593,018 rows | 4,862 products x 739 days | zero-share 86% | median price GBP 2.46
2 · The long tail¶
The full catalogue is a different animal from the best-sellers: most products sell rarely. This is the regime a real planner faces — a handful of stars and thousands of slow movers.
In [3]:
active = c.groupby("StockCode")["date"].nunique(); total = c.groupby("StockCode")["Quantity"].sum()
print(f"active-days per product: median {active.median():.0f} of {len(all_days)} | p25 {active.quantile(.25):.0f} | p75 {active.quantile(.75):.0f}")
print(f"sell on <30 days: {(active<30).mean():.0%} | top-200 hold {total.sort_values(ascending=False).head(200).sum()/total.sum():.0%} of all demand")
fig, ax = plt.subplots(1, 2, figsize=(13, 4))
ax[0].hist(active, bins=60, color=BLUE); ax[0].axvline(active.median(), color=RED, ls="--", label=f"median {active.median():.0f}")
ax[0].set_xlabel("active (nonzero) days"); ax[0].set_ylabel("# products"); ax[0].set_title("Most products sell on few days"); ax[0].legend(fontsize=8)
sv = total.sort_values(ascending=False).values; cumsh = np.cumsum(sv)/sv.sum()
ax[1].plot(np.arange(1, len(cumsh)+1), cumsh, color=GREEN); ax[1].axvline(200, color=RED, ls="--", label="top 200")
ax[1].axhline(cumsh[199], color=GREY, ls=":"); ax[1].set_xscale("log")
ax[1].set_xlabel("products (ranked by volume, log)"); ax[1].set_ylabel("cumulative share of demand")
ax[1].set_title("A few products carry most demand"); ax[1].legend(fontsize=8)
fig.tight_layout(); plt.show()
active-days per product: median 62 of 739 | p25 20 | p75 145 sell on <30 days: 33% | top-200 hold 43% of all demand
3 · Demand character — now genuinely intermittent¶
The Syntetos–Boylan lens (see the top-200 notebook §3 for the Croston model and the 1.32 / 0.49 cut-offs). Where the best-sellers were almost all lumpy (sell often, variable size), the full catalogue brings in the intermittent slow movers (sparse but steady size) — the classic intermittent-demand regime.
In [4]:
def sb_stats(df):
def one(g):
y = g["demand"].to_numpy(); nz = y[y>0]
return pd.Series({"mean": y.mean(), "zero_share": (y==0).mean(),
"adi": len(y)/max(len(nz),1), "cv2": (nz.std()/nz.mean())**2 if len(nz)>1 else 0.0})
return df.groupby("StockCode").apply(one, include_groups=False)
def sb_class(a,c2):
if a<1.32 and c2<0.49: return "smooth"
if a>=1.32 and c2<0.49: return "intermittent"
if a<1.32 and c2>=0.49: return "erratic"
return "lumpy"
ps = sb_stats(panel); ps["class"] = [sb_class(a,c2) for a,c2 in zip(ps.adi, ps.cv2)]
print((ps["class"].value_counts(normalize=True)*100).round(1).to_string())
cls_col = {"smooth":GREEN,"intermittent":BLUE,"erratic":ORANGE,"lumpy":RED}
fig, ax = plt.subplots(1, 2, figsize=(13, 4.5))
for cl,col in cls_col.items():
m = ps["class"]==cl
ax[0].scatter(ps.adi[m], ps.cv2[m], s=6, color=col, alpha=.4, label=cl)
ax[0].axvline(1.32, color=GREY, ls="--"); ax[0].axhline(0.49, color=GREY, ls="--")
ax[0].set_xlim(1, min(ps.adi.quantile(.97), 30)); ax[0].set_ylim(0, min(ps.cv2.quantile(.97), 4))
ax[0].set_xlabel("ADI"); ax[0].set_ylabel("CV² of nonzero demand"); ax[0].set_title("Syntetos–Boylan (all 4,862 products)"); ax[0].legend(fontsize=8)
vc = ps["class"].value_counts()
ax[1].bar(vc.index, vc.values, color=[cls_col[k] for k in vc.index]); ax[1].set_title("Products by class"); ax[1].set_ylabel("# products")
fig.tight_layout(); plt.show()
class lumpy 82.6 intermittent 17.2 erratic 0.2
Intermittent demand — sporadic but consistent size — now makes up a real share (the slow-mover tail), alongside the lumpy best-sellers. Both are regimes where a point forecast collapses toward zero and a distribution is required. Zero-share across the panel is ~86%, so the probability atom at zero will dominate the low quantiles even more than before.
4 · Features — the same leak-safe matrix, vectorized for scale¶
Identical construction to the top-200 notebook §4 (the avail = demand.shift(H) primitive keeps every
history feature leak-safe), but the days_since_last_sale computation is vectorized so it runs on 3.6M
rows in seconds rather than looping per product.
In [5]:
HORIZON, LAGS, WINDOWS, SERIES = 7, [7,14,28], [7,28], "StockCode"
def build_features(panel, ps, H, lags, windows):
d = panel.sort_values([SERIES,"date"]).reset_index(drop=True).copy()
dt = d["date"].dt
d["dow"]=dt.dayofweek; d["is_saturday"]=(d["dow"]==5).astype("int8"); d["month"]=dt.month; d["day"]=dt.day
d["weekofyear"]=dt.isocalendar().week.astype("int32"); d["dayofyear"]=dt.dayofyear
ang=2*np.pi*d["dayofyear"]/365.25; d["sin_year"]=np.sin(ang); d["cos_year"]=np.cos(ang)
dxm=(d["dayofyear"]-359).abs(); d["days_to_xmas"]=np.minimum(dxm,365-dxm)
d["log_price"]=np.log(d["price"].clip(lower=1e-3))
feats=["dow","is_saturday","month","day","weekofyear","dayofyear","sin_year","cos_year","days_to_xmas","price","log_price"]
avail=d.groupby(SERIES, observed=True)["demand"].shift(H); gav=avail.groupby(d[SERIES], observed=True)
d["lag_0"]=avail; feats.append("lag_0")
for k in lags: d[f"lag_{k}"]=gav.shift(k); feats.append(f"lag_{k}")
for w in windows:
d[f"roll_mean_{w}"]=gav.transform(lambda s: s.rolling(w,min_periods=1).mean())
d[f"roll_std_{w}"]=gav.transform(lambda s: s.rolling(w,min_periods=2).std())
d[f"roll_zero_{w}"]=gav.transform(lambda s: s.eq(0).rolling(w,min_periods=1).mean())
feats+=[f"roll_mean_{w}",f"roll_std_{w}",f"roll_zero_{w}"]
gc=d.groupby(SERIES, observed=True).cumcount(); marker=np.where((avail>0).to_numpy(), gc.to_numpy(), np.nan)
last=pd.Series(marker, index=d.index).groupby(d[SERIES], observed=True).ffill()
d["days_since_last_sale"]=gc.to_numpy()-last.to_numpy(); feats.append("days_since_last_sale")
d["sb_class"]=d[SERIES].map(ps["class"]).astype("category"); d[SERIES]=d[SERIES].astype("category")
feats+=[SERIES,"sb_class"]
return d, feats, [SERIES,"sb_class"]
t=time.time(); feat_df, FEATURES, CAT = build_features(panel, ps, HORIZON, LAGS, WINDOWS)
print(f"feature matrix {len(feat_df):,} rows x {len(FEATURES)} features built in {time.time()-t:.1f}s")
feature matrix 3,593,018 rows x 24 features built in 5.6s
5 · Distributional forecast (quantile LightGBM)¶
Same quantile-boosting model as the top-200 notebook §5 (one LightGBM per quantile, pinball loss, non-crossing). It scales to 3.6M rows in about a minute. We check calibration — with 86% zeros, the low-quantile atom is even more pronounced.
In [6]:
import lightgbm as lgb
QUANTILES=[0.05,0.1,0.25,0.5,0.75,0.9,0.95,0.975,0.99]
PARAMS=dict(objective="quantile",n_estimators=300,learning_rate=0.05,num_leaves=63,min_child_samples=50,
subsample=0.8,subsample_freq=1,colsample_bytree=0.8,verbosity=-1,n_jobs=-1)
def pinball(y,q,tau): e=np.asarray(y,float)-np.asarray(q,float); return float(np.mean(np.maximum(tau*e,(tau-1)*e)))
class QuantileGBM:
def __init__(s,qs,p): s.q=sorted(qs); s.p=p; s.m={}
def fit(s,X,y,Xv,yv):
for t in s.q:
mm=lgb.LGBMRegressor(alpha=t,random_state=7,**s.p)
mm.fit(X,y,eval_set=[(Xv,yv)],eval_metric="quantile",callbacks=[lgb.early_stopping(40,verbose=False),lgb.log_evaluation(0)]); s.m[t]=mm
return s
def predict(s,X):
P=np.clip(np.column_stack([s.m[t].predict(X) for t in s.q]),0,None); P.sort(axis=1)
return pd.DataFrame(P,columns=s.q,index=X.index)
cut=feat_df["date"].max()-pd.Timedelta(days=28); tr=feat_df[feat_df.date<=cut]; va=feat_df[feat_df.date>cut]
t=time.time(); model=QuantileGBM(QUANTILES,PARAMS).fit(tr[FEATURES],tr["demand"],va[FEATURES],va["demand"])
zshare=float((va.demand==0).mean())
val=pd.DataFrame([dict(quantile=q,coverage=float((va.demand.to_numpy()<=model.predict(va[FEATURES])[q].to_numpy()).mean())) for q in QUANTILES])
print(f"trained 9 quantile models on {len(tr):,} rows in {time.time()-t:.0f}s | validation zero-share {zshare:.0%}")
fig,ax=plt.subplots(figsize=(6,5))
ax.plot([0,1],[0,1],"--",color=GREY,label="perfect"); ax.plot(val["quantile"],val["coverage"],"o-",color=BLUE,label="empirical")
ax.axhline(zshare,color=RED,ls=":",label=f"zero-share {zshare:.2f}")
ax.set_xlabel("target quantile τ"); ax.set_ylabel("coverage"); ax.set_title("Calibration — full catalogue"); ax.legend(fontsize=8)
fig.tight_layout(); plt.show()
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)
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)
trained 9 quantile models on 3,456,882 rows in 74s | validation zero-share 77%
As expected, the flat low-quantile arm now sits near 0.86 (the panel zero-share): for products that are zero ~86% of days, every quantile up to ~0.86 is genuinely 0. The upper tail — where the newsvendor reads — is what must be calibrated, and it is.
6 · The decision — newsvendor with per-product economics¶
Same economics as the top-200 notebook §6 ($C_u=0.5\times$ price, $C_o=£0.10$ flat → a per-product critical fractile that rises with price). On the long tail, most orders are small — often 0 or 1 — which is exactly the point: the model decides which slow movers are worth stocking at all.
In [7]:
MARGIN_RATE, HOLDING = 0.5, 0.10
price=panel.groupby("StockCode")["price"].first(); cu_prod=MARGIN_RATE*price; fractile=cu_prod/(cu_prod+HOLDING)
def order_at_fractile(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))])
snap=feat_df[feat_df["date"]==feat_df["date"].max()].copy()
Psnap=model.predict(snap[FEATURES]); Psnap.index=snap["StockCode"].to_numpy()
fr=fractile.reindex(Psnap.index).to_numpy(); order=pd.Series(order_at_fractile(Psnap,fr,QUANTILES),index=Psnap.index)
print(f"snapshot {snap['date'].iloc[0].date()} | {len(order):,} products | total order {order.sum():,.0f} units")
print(f"products ordered 0 units: {(order<0.5).mean():.0%} | ordered >=1: {(order>=0.5).sum():,}")
fig,ax=plt.subplots(figsize=(7,4))
ax.hist(order.clip(upper=order.quantile(.99)),bins=60,color=GREEN); ax.set_yscale("log")
ax.set_xlabel("recommended order (units, snapshot day)"); ax.set_ylabel("# products (log)")
ax.set_title("Most of the long tail gets a small order — a few get large ones"); fig.tight_layout(); plt.show()
snapshot 2011-12-09 | 4,862 products | total order 51,408 units products ordered 0 units: 43% | ordered >=1: 2,782
7 · Backtest — the value of modeling uncertainty on the full catalogue¶
Rolling-origin, three policies (point-mean / median / per-product newsvendor), realized £ cost with per-product economics — as in the top-200 notebook §7. The long tail is where point forecasts fail hardest, so we expect the newsvendor advantage to be at least as large.
In [8]:
BT_Q=[0.25,0.5,0.7,0.8,0.9,0.95,0.99]; N_FOLDS, FOLD_H=3,28
def fold_windows(mx,n,h): return list(reversed([(mx-pd.Timedelta(days=k*h)-pd.Timedelta(days=h-1), mx-pd.Timedelta(days=k*h)) for k in range(n)]))
def run_bt():
parts=[]
for fold,(ts,te) in enumerate(fold_windows(feat_df.date.max(),N_FOLDS,FOLD_H)):
trn=feat_df[feat_df.date<ts]; tst=feat_df[(feat_df.date>=ts)&(feat_df.date<=te)].copy()
cv=trn.date.max()-pd.Timedelta(days=FOLD_H); f2,v2=trn[trn.date<=cv],trn[trn.date>cv]
qm=QuantileGBM(BT_Q,PARAMS).fit(f2[FEATURES],f2["demand"],v2[FEATURES],v2["demand"])
mp=dict(PARAMS); mp["objective"]="regression"; mm=lgb.LGBMRegressor(random_state=7,**mp)
mm.fit(f2[FEATURES],f2["demand"],eval_set=[(v2[FEATURES],v2["demand"])],eval_metric="l2",callbacks=[lgb.early_stopping(40,verbose=False),lgb.log_evaluation(0)])
P=qm.predict(tst[FEATURES]); frr=fractile.reindex(tst["StockCode"]).to_numpy()
tst["Q_newsvendor"]=order_at_fractile(P,frr,BT_Q); tst["Q_median"]=P[0.5].to_numpy()
tst["Q_point"]=np.clip(mm.predict(tst[FEATURES]),0,None); tst["cu"]=cu_prod.reindex(tst["StockCode"]).to_numpy(); tst["fold"]=fold
parts.append(tst[["fold","date","StockCode","demand","cu","Q_point","Q_median","Q_newsvendor"]])
return pd.concat(parts,ignore_index=True)
t=time.time(); bt=run_bt(); D=bt.demand.to_numpy(float); CU=bt.cu.to_numpy(float)
POL={"point (mean)":"Q_point","median (q0.5)":"Q_median","newsvendor":"Q_newsvendor"}
for n,qc in POL.items(): bt[qc.replace("Q_","cost_")]=CU*np.maximum(D-bt[qc].to_numpy(float),0)+HOLDING*np.maximum(bt[qc].to_numpy(float)-D,0)
tot={n:bt[qc.replace("Q_","cost_")].sum() for n,qc in POL.items()}
voi=(tot["point (mean)"]-tot["newsvendor"])/tot["point (mean)"]
print(f"backtest {N_FOLDS}x{FOLD_H}d on {len(bt):,} product-days in {time.time()-t:.0f}s")
print(f"realized cost: point GBP {tot['point (mean)']:,.0f} | median GBP {tot['median (q0.5)']:,.0f} | newsvendor GBP {tot['newsvendor']:,.0f}")
print(f"VALUE OF MODELING UNCERTAINTY: newsvendor cuts cost {voi:+.1%} vs point-forecast policy")
COL={"point (mean)":ORANGE,"median (q0.5)":PURP,"newsvendor":GREEN}
daily=bt.groupby("date")[[qc.replace("Q_","cost_") for qc in POL.values()]].sum().cumsum()
fig,ax=plt.subplots(figsize=(12,5))
for n,qc in POL.items(): ax.plot(daily.index,daily[qc.replace("Q_","cost_")],lw=2,color=COL[n],label=n)
ax.set_ylabel("cumulative realized cost (GBP)"); ax.set_title("Realized cost over the backtest — full catalogue"); ax.legend(fontsize=9)
fig.tight_layout(); plt.show()
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)
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)
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)
backtest 3x28d on 408,408 product-days in 190s realized cost: point GBP 1,089,435 | median GBP 1,246,814 | newsvendor GBP 822,435 VALUE OF MODELING UNCERTAINTY: newsvendor cuts cost +24.5% vs point-forecast policy
8 · Wrap-up — the tail is where it pays¶
Run on all 4,862 products, the pipeline is unchanged but the stakes shift: the catalogue is ~86% zeros, dominated by intermittent and lumpy slow movers, and that is exactly where a point forecast rounds to "stock nothing." The distribution-driven newsvendor policy again cuts realized cost against an equally-good point forecaster — the dollar value of modeling uncertainty, now measured across the whole catalogue.
Two honest limits this scale exposes:
- Ultra-slow movers (sold on a handful of days) carry almost no signal — collectively they matter for capital tied up, individually the decision is near-trivial (stock 0 or 1).
- Cold-start / sparse products are where a pooling model — a Bayesian hierarchical hurdle with product random effects — should beat a global GBM, because priors and partial pooling substitute for missing history. That is the subject of the third notebook in this folder, which races the two on a stratified sample by data-richness.
Together the three notebooks make the full argument: the method is sound on the best-sellers (notebook 1), it scales to the whole catalogue (this notebook), and where data is scarce, a Bayesian pooling model is the competitor worth taking seriously (notebook 3).