Staffing under uncertainty — forecasting meets a two-stage stochastic program¶
Machine Learning in Operations Research · decision under uncertainty, done properly¶
The predict-then-optimize notebooks solved a one-stage decision (the newsvendor: order once, demand realizes, done). This one takes the next rung: a two-stage stochastic program — the OR workhorse for decisions committed before uncertainty resolves, with recourse afterwards — and it does the full arc, classical → machine learning → Bayesian → OR integration.
The problem. An emergency department must provision capacity (nurses / bed-hours) for the coming week before it knows how many patients will arrive. Under-provision → expensive overtime / agency cover; over-provision → wasted committed cost. A fixed weekly capacity budget must be allocated across the seven days, so the days compete and the allocation should hedge toward days that are both busy and uncertain.
The plan of attack.
- Establish the demand and its real seasonality (§1–2).
- Forecast it with the methods from the paper this data came from — ARIMA, Holt-Winters, TBATS, a neural network — plus a Bayesian model, and compare their accuracy (§3).
- Argue that a point forecast is the wrong object for an asymmetric decision, and turn the Bayesian model's posterior predictive into a scenario set (§4, §6).
- Formulate and solve the two-stage stochastic program, and measure — in $ — the Value of the Stochastic Solution (VSS) and the Expected Value of Perfect Information (EVPI) (§5, §7–8).
The data is real: hourly ED arrivals at UnityPoint Health (Iowa), 2014–2017, ~100 patients/day (Dryad); the source study benchmarked ARIMA, Holt-Winters, TBATS and neural networks for the forecasting task.
import os
os.environ.setdefault("KMP_DUPLICATE_LIB_OK", "TRUE")
import numpy as np, pandas as pd, matplotlib.pyplot as plt, warnings
from pathlib import Path
warnings.filterwarnings("ignore")
BLUE, RED, GREEN, ORANGE, GREY, PURP, TEAL = "#2b6cb0","#c53030","#2f855a","#dd6b20","#a0aec0","#6b46c1","#0891b2"
here = Path.cwd(); DATA = next(p for p in [here/"data", here.parent/"data"] if (p/"ed_hourly.parquet").exists())
h = pd.read_parquet(DATA/"ed_hourly.parquet")["arrivals"]; h.index = pd.to_datetime(h.index)
daily = h.resample("D").sum(); daily = daily[daily >= 30] # drop partial edge days (mean ~100)
print(f"hourly obs {len(h):,} | daily obs {len(daily):,} | {daily.index.min().date()} -> {daily.index.max().date()}")
print(f"daily arrivals: mean {daily.mean():.1f}, std {daily.std():.1f}, min {daily.min():.0f}, max {daily.max():.0f}")
hourly obs 33,984 | daily obs 1,412 | 2014-01-01 -> 2017-12-30 daily arrivals: mean 100.6, std 15.9, min 54, max 153
1 · The data — real emergency-department arrivals¶
One row per calendar day: the number of patients arriving at the ED. Unlike a synthetic demo, this is a real operational series, so it carries the structure a staffing plan must respect — and the irreducible noise it must hedge against. We drop a few partial edge days (counts far below the ~100/day norm).
2 · Demand character — three layers of real seasonality¶
A staffing plan lives or dies on seasonality. There are three layers here: a time-of-day profile (who arrives when), a day-of-week profile (which days are busy), and a slow multi-year trend. A good forecast captures these on average; a good stochastic plan additionally hedges the spread around them.
fig, ax = plt.subplots(2, 2, figsize=(13, 8))
ax[0,0].plot(daily.index, daily.values, color=GREY, lw=.4); ax[0,0].plot(daily.index, daily.rolling(28).mean(), color=RED, lw=2)
ax[0,0].set_title("Daily ED arrivals (28-day mean): trend"); ax[0,0].set_ylabel("patients/day")
byh = h.groupby(h.index.hour).mean()
ax[0,1].plot(byh.index, byh.values, "o-", color=BLUE); ax[0,1].set_title("Time-of-day profile (real diurnal shape)"); ax[0,1].set_xlabel("hour"); ax[0,1].set_ylabel("arrivals/hour")
names = ["Mon","Tue","Wed","Thu","Fri","Sat","Sun"]
bydow = daily.groupby(daily.index.dayofweek).mean()
ax[1,0].bar(range(7), bydow.values, color=GREEN); ax[1,0].set_xticks(range(7)); ax[1,0].set_xticklabels(names)
ax[1,0].set_title("Day-of-week profile (Monday peak)"); ax[1,0].set_ylabel("mean patients/day"); ax[1,0].set_ylim(bydow.min()-5, bydow.max()+3)
bymo = daily.groupby(daily.index.month).mean()
ax[1,1].plot(bymo.index, bymo.values, "o-", color=PURP); ax[1,1].set_xticks(range(1,13)); ax[1,1].set_title("Month-of-year"); ax[1,1].set_xlabel("month")
fig.suptitle("Real ED arrival seasonality — UnityPoint Health 2014–2017", fontsize=13); fig.tight_layout(); plt.show()
print(f"day-of-week spread: {bydow.min():.0f} (Thu) to {bydow.max():.0f} (Mon) patients/day | diurnal peak ~{byh.max():.1f}/hr vs trough ~{byh.min():.1f}/hr")
day-of-week spread: 96 (Thu) to 107 (Mon) patients/day | diurnal peak ~6.0/hr vs trough ~1.6/hr
3 · Forecasting the arrivals — a model bake-off¶
The study that published this data benchmarked four forecasters. We reproduce that comparison and add a Bayesian model, because — as §4 will argue — the decision needs a distribution, which only the Bayesian model yields natively. All are fit on the first 75% of days and forecast the held-out 25%.
- ARIMA / SARIMA — AutoRegressive Integrated Moving Average. Models the series through its own autocorrelation and moving-average errors; the seasonal form $\text{SARIMA}(p,d,q)(P,D,Q)_7$ adds a weekly cycle. The classical time-series default.
- Holt-Winters (ETS) — exponential smoothing with additive trend and additive weekly seasonal components, each updated by a smoothing recursion. Simple, robust, hard to beat on clean seasonal data.
- TBATS — Trigonometric seasonality, Box-Cox, ARMA errors, Trend, Seasonal. Represents seasonality with Fourier terms and layers an ARMA error model on top; built for complex/multiple seasonalities.
- Neural network — a small multilayer perceptron mapping calendar features (day-of-week, month, trend, annual harmonics) to arrivals; the ML/DL representative.
- Bayesian NegBin GLM — a Negative-Binomial regression on the same calendar features, fit by MCMC. Alone among these it returns a full posterior predictive distribution, not just a point.
ntr = int(len(daily)*0.75); tr, te = daily.iloc[:ntr], daily.iloc[ntr:]
y_te = te.values; fc = {}
# --- calendar feature frame (for NN + Bayesian) ---
def feats(idx, pos):
doy = idx.dayofyear.values
return pd.DataFrame({"dow": idx.dayofweek.values, "month": idx.month.values-1,
"t": (pos - len(daily)/2)/len(daily),
"sin": np.sin(2*np.pi*doy/365.25), "cos": np.cos(2*np.pi*doy/365.25)}, index=idx)
F = feats(daily.index, np.arange(len(daily))); Ftr, Fte = F.iloc[:ntr], F.iloc[ntr:]
# --- ARIMA (SARIMA weekly) ---
from statsmodels.tsa.statespace.sarimax import SARIMAX
fc["ARIMA"] = SARIMAX(tr.values, order=(2,0,1), seasonal_order=(1,0,1,7),
enforce_stationarity=False, enforce_invertibility=False).fit(disp=False).forecast(len(te))
# --- Holt-Winters ---
from statsmodels.tsa.holtwinters import ExponentialSmoothing
fc["Holt-Winters"] = ExponentialSmoothing(tr.values, trend="add", seasonal="add", seasonal_periods=7).fit().forecast(len(te))
# --- TBATS (one fit; ~50s) ---
from tbats import TBATS
fc["TBATS"] = TBATS(seasonal_periods=[7], use_box_cox=False, use_arma_errors=True, n_jobs=1).fit(tr.values).forecast(len(te))
print("classical models done")
C:\Users\user\anaconda3\envs\pymc-env\Lib\site-packages\statsmodels\base\model.py:607: ConvergenceWarning: Maximum Likelihood optimization failed to converge. Check mle_retvals
warnings.warn("Maximum Likelihood optimization failed to "
classical models done
# --- Neural network (MLP on calendar features) ---
import torch, torch.nn as nn
def Xmat(Fsub):
dow = np.eye(7)[Fsub["dow"].values]; mon = np.eye(12)[Fsub["month"].values]
cont = Fsub[["t","sin","cos"]].values
return np.hstack([dow, mon, cont]).astype(np.float32)
Xtr, Xte = Xmat(Ftr), Xmat(Fte)
ymu, ysd = tr.mean(), tr.std()
net = nn.Sequential(nn.Linear(Xtr.shape[1], 32), nn.ReLU(), nn.Linear(32, 1))
opt = torch.optim.Adam(net.parameters(), lr=1e-2); Xt = torch.tensor(Xtr); yt = torch.tensor(((tr.values-ymu)/ysd).astype(np.float32))
for ep in range(400):
opt.zero_grad(); loss = ((net(Xt).squeeze()-yt)**2).mean(); loss.backward(); opt.step()
with torch.no_grad(): fc["Neural net"] = net(torch.tensor(Xte)).squeeze().numpy()*ysd + ymu
# --- Bayesian NegBin GLM (fit once; reused for scenarios later) ---
import jax, jax.numpy as jnp, numpyro, numpyro.distributions as dist
from numpyro.infer import MCMC, NUTS
def model(dow, month, t, sin, cos, y=None):
a = numpyro.sample("a", dist.Normal(np.log(daily.mean()), 1.0))
b_dow = numpyro.sample("b_dow", dist.Normal(0,0.5).expand([7]))
b_mon = numpyro.sample("b_mon", dist.Normal(0,0.5).expand([12]))
b_t = numpyro.sample("b_t", dist.Normal(0,0.5)); b_s = numpyro.sample("b_s", dist.Normal(0,0.5).expand([2]))
conc = numpyro.sample("conc", dist.HalfNormal(50.))
mu = jnp.exp(a + b_dow[dow] + b_mon[month] + b_t*t + b_s[0]*sin + b_s[1]*cos)
numpyro.sample("obs", dist.NegativeBinomial2(mu, conc), obs=y)
args_tr = (jnp.array(Ftr.dow.values), jnp.array(Ftr.month.values), jnp.array(Ftr.t.values), jnp.array(Ftr.sin.values), jnp.array(Ftr.cos.values))
mcmc = MCMC(NUTS(model), num_warmup=500, num_samples=600, num_chains=1, progress_bar=False)
mcmc.run(jax.random.PRNGKey(0), *args_tr, y=jnp.array(tr.values, float)); post = mcmc.get_samples()
def bayes_mu(Fsub):
return np.exp(post["a"].mean() + post["b_dow"].mean(0)[Fsub.dow.values] + post["b_mon"].mean(0)[Fsub.month.values]
+ post["b_t"].mean()*Fsub.t.values + post["b_s"].mean(0)[0]*Fsub.sin.values + post["b_s"].mean(0)[1]*Fsub.cos.values)
fc["Bayesian"] = bayes_mu(Fte)
print("neural net + Bayesian done")
neural net + Bayesian done
def metrics(a, f):
e = a-f; return dict(MAE=np.mean(np.abs(e)), RMSE=np.sqrt(np.mean(e**2)), MAPE=np.mean(np.abs(e)/a)*100)
tbl = pd.DataFrame({k: metrics(y_te, v) for k, v in fc.items()}).T.sort_values("MAE")
print("held-out forecast accuracy (25% of days):"); print(tbl.round(2).to_string())
fig, ax = plt.subplots(figsize=(13, 4.5))
sl = slice(0, 70)
ax.plot(te.index[sl], y_te[sl], color="black", lw=1.6, label="actual", zorder=5)
for (k, v), c in zip(fc.items(), [BLUE, GREEN, ORANGE, PURP, RED]):
ax.plot(te.index[sl], v[sl], lw=1.2, color=c, alpha=.8, label=k)
ax.set_title("Forecasts vs actual (first 70 held-out days)"); ax.set_ylabel("patients/day"); ax.legend(fontsize=8, ncol=3); fig.tight_layout(); plt.show()
held-out forecast accuracy (25% of days):
MAE RMSE MAPE
Bayesian 10.14 12.69 8.87
Holt-Winters 10.45 12.85 9.44
TBATS 10.48 12.94 9.28
ARIMA 10.80 13.16 9.48
Neural net 12.11 15.04 10.88
Reading the bake-off. All five sit within a narrow band — strong, stable seasonality is easy to capture, so the choice of forecaster barely moves point accuracy (the source paper found the same: the methods cluster). That is the important lesson going in: on a clean seasonal series, a better point forecast is not where the value is. The value is in what you do with the forecast — and in carrying its uncertainty into that decision, which is where the Bayesian model and the optimizer come in.
4 · Why a point forecast is the wrong object¶
Suppose the best model says Monday will see 110 patients. Staffing for exactly 110 is only optimal if being one patient short costs the same as being one over. It does not: a shortfall pulls in overtime / agency cover at a premium, while slack is merely the (smaller) committed cost. With this asymmetry the right capacity is not the mean — it is a higher quantile of the demand distribution. So the decision needs the whole distribution, and a plan built on the point forecast will systematically under-provision.
This is exactly the newsvendor logic of the predict-then-optimize notebooks — and it generalizes here to a two-stage problem, because capacity is committed across the week under a shared budget before any of the seven days' demand is known.
5 · The two-stage stochastic program¶
Let $C_d\ge 0$ be the capacity provisioned for day $d$ (patients handled at regular cost) — the first-stage decision, made before arrivals are known, under a fixed weekly capacity budget $B$:
$$\sum_{d=1}^{7} C_d \le B,\qquad C_d\ge 0.$$
Then weekly arrivals $A=(A_1,\dots,A_7)$ realize and the second-stage recourse covers the shortfall $u_d=(A_d-C_d)^+$ with overtime at a higher unit cost. With regular unit cost $c_r$ and overtime $c_o>c_r$:
$$\min_{C\ge 0,\ \sum_d C_d\le B}\ \ c_r\sum_d C_d\ +\ \mathbb{E}_A\!\Big[\ \min_{u\ge 0,\ u_d\ge A_d-C_d}\ c_o\sum_d u_d\ \Big].$$
The expectation is over a scenario set $\{A^{(s)}\}_{s=1}^S$, which turns the program into one linear program. Two ways to build that scenario set define the contrast we measure:
- Deterministic (EEV) — replace $A$ by its mean $\bar A$ (the point forecast). The naive plan.
- Stochastic (RP) — use $S$ scenarios that span the spread; from a Bayesian posterior predictive they carry parameter and sampling uncertainty.
We solve with PuLP/CBC. Cost units are illustrative ($c_r=1$, overtime $c_o=3\times$).
import pulp
def solve_sp(scen, budget, c_r=1.0, c_o=3.0):
S, D = scen.shape
prob = pulp.LpProblem("staffing", pulp.LpMinimize)
C = {d: pulp.LpVariable(f"C{d}", lowBound=0) for d in range(D)}
u = {(s,d): pulp.LpVariable(f"u{s}_{d}", lowBound=0) for s in range(S) for d in range(D)}
prob += c_r*pulp.lpSum(C[d] for d in range(D)) + (c_o/S)*pulp.lpSum(u[s,d] for s in range(S) for d in range(D))
for s in range(S):
for d in range(D):
prob += u[s,d] >= float(scen[s,d]) - C[d]
prob += pulp.lpSum(C[d] for d in range(D)) <= budget
prob.solve(pulp.PULP_CBC_CMD(msg=0))
return np.array([C[d].value() or 0.0 for d in range(D)])
def realized_cost(C, actual, c_r=1.0, c_o=3.0):
return c_r*np.sum(C) + c_o*np.sum(np.maximum(actual - C, 0))
print("SP solver ready (PuLP/CBC)")
SP solver ready (PuLP/CBC)
6 · Scenarios from the Bayesian posterior predictive¶
The Bayesian NegBin model already fitted gives, for any future day, a full predictive distribution. To build a scenario for an upcoming week we draw one posterior sample of the parameters, compute each day's mean, and draw a Negative-Binomial count — repeating $S$ times. The result is a cloud of plausible weeks that respects day-of-week, seasonal, trend and the overdispersed noise. A posterior-predictive check confirms the model reproduces the demand's spread.
def ppred(Fsub, n=200, key=1):
idx = np.random.default_rng(key).choice(len(post["a"]), n, replace=False)
mu = np.exp(post["a"][idx][:,None] + post["b_dow"][idx][:,Fsub.dow.values] + post["b_mon"][idx][:,Fsub.month.values]
+ post["b_t"][idx][:,None]*Fsub.t.values[None,:] + post["b_s"][idx][:,0][:,None]*Fsub.sin.values[None,:]
+ post["b_s"][idx][:,1][:,None]*Fsub.cos.values[None,:])
r = post["conc"][idx][:,None]; p = r/(r+mu)
return np.random.default_rng(key).negative_binomial(r, p)
pp = ppred(Ftr, n=400)
fig, ax = plt.subplots(1, 2, figsize=(13, 4))
ax[0].hist(tr.values, bins=40, density=True, alpha=.6, color=GREY, label="actual (train)")
ax[0].hist(pp.reshape(-1)[np.random.default_rng(0).integers(0, pp.size, 20000)], bins=40, density=True, histtype="step", color=RED, lw=2, label="posterior predictive")
ax[0].set_title("Posterior predictive check"); ax[0].set_xlabel("daily arrivals"); ax[0].legend(fontsize=8)
ax[1].bar(range(7), np.exp(post["b_dow"].mean(0))-1, color=BLUE); ax[1].set_xticks(range(7)); ax[1].set_xticklabels(names)
ax[1].axhline(0, color=GREY, lw=.8); ax[1].set_title("Day-of-week effect (posterior mean, ×average)"); ax[1].set_ylabel("relative to average")
fig.tight_layout(); plt.show()
Reading the two panels¶
Left — the posterior predictive check. This is the model auditing itself as a scenario generator. The grey histogram is the actual distribution of daily arrivals in the training period; the red outline is the distribution the fitted model simulates (draw parameters from the posterior, then draw a Negative-Binomial count for each). They overlay closely — same centre (~100/day), same spread, same mild right-skew — which is exactly what we need: the model reproduces not just average demand but its variability and shape, so the weeks it hands the optimizer look like real weeks. (A model that matched the mean but under-dispersed would generate over-confident scenarios and an under-hedged plan — the check is there to rule that out.)
Right — the day-of-week effect. Each bar is the posterior-mean multiplier for a weekday relative to the overall average (0 = average; +0.07 ≈ 7% busier than a typical day). Monday sits highest — the classic ED Monday peak — with a midweek dip (Thursday lowest) and a weekend rebound, mirroring the raw profile in §2 (Monday ~107 vs Thursday ~95 patients/day). This is the systematic structure every scenario inherits: when we simulate a week, its Monday is centred high and its Thursday low, so the stochastic plan provisions more capacity for Mondays. Crucially, around each of these day-specific centres the posterior predictive still spreads the count noise — and that residual spread, not the seasonal mean, is the uncertainty the two-stage program exists to hedge.
7 · Solve — deterministic vs stochastic, on one week¶
For a held-out week we draw the scenario set and solve the SP both ways. The stochastic plan provisions above the mean forecast and tilts capacity toward the busier, more variable days — the hedge a point forecast cannot express, all within the same tight budget.
te_feat = feats(te.index, np.arange(ntr, len(daily)))
weeks_idx = [g for _, g in te.groupby(te.index.isocalendar().week.astype(str)+"-"+te.index.year.astype(str)) if len(g)==7]
BUDGET = 0.95 * 7 * daily.mean()
exw = weeks_idx[3]; exF = te_feat.loc[exw.index]
scen = ppred(exF, n=200); mean_fc = scen.mean(0)
C_det = solve_sp(mean_fc[None,:], BUDGET); C_sto = solve_sp(scen, BUDGET)
fig, ax = plt.subplots(figsize=(11, 4.5)); x = np.arange(7); w = 0.28
ax.bar(x-w, mean_fc, w, color=GREY, label="mean forecast (demand)")
ax.bar(x, C_det, w, color=ORANGE, label="deterministic plan (EEV)")
ax.bar(x+w, C_sto, w, color=GREEN, label="stochastic plan (RP)")
ax.plot(x, exw.values, "o", color=RED, ms=9, label="actual arrivals")
ax.set_xticks(x); ax.set_xticklabels([f"{n}\n{d.strftime('%m/%d')}" for n, d in zip(names, exw.index)])
ax.set_ylabel("patients"); ax.set_title(f"Week of {exw.index.min().date()} — capacity plans vs demand"); ax.legend(fontsize=8)
fig.tight_layout(); plt.show()
print(f"budget B={BUDGET:.0f} | deterministic realized ${realized_cost(C_det, exw.values):.0f} | stochastic realized ${realized_cost(C_sto, exw.values):.0f}")
budget B=669 | deterministic realized $1122 | stochastic realized $1062
8 · Backtest — VSS and EVPI over every held-out week¶
Rolling the plan over all complete held-out weeks and accumulating realized cost gives the two numbers OR practitioners look for:
- RP (stochastic, Bayesian scenarios), EEV (deterministic mean plan), WS (wait-and-see: the impossible plan that knows demand in advance — a lower bound).
- VSS $=$ EEV $-$ RP — what the stochastic solution saves over the naive mean plan (the value of modelling uncertainty, in $).
- EVPI $=$ RP $-$ WS — what perfect information would additionally be worth.
rows = []
for wk in weeks_idx:
wF = te_feat.loc[wk.index]; sc = ppred(wF, n=200); a = wk.values
rows.append(dict(week=wk.index.min().date(),
RP=realized_cost(solve_sp(sc, BUDGET), a),
EEV=realized_cost(solve_sp(sc.mean(0)[None,:], BUDGET), a),
WS=realized_cost(solve_sp(a[None,:], BUDGET), a)))
R = pd.DataFrame(rows); tot = R[["EEV","RP","WS"]].sum()
VSS = tot["EEV"]-tot["RP"]; EVPI = tot["RP"]-tot["WS"]
wins = int((R["RP"] <= R["EEV"]).sum())
print(f"held-out weeks: {len(R)} | RP <= EEV in {wins}/{len(R)} weeks")
print(f"total realized cost EEV ${tot['EEV']:,.0f} | RP ${tot['RP']:,.0f} | WS ${tot['WS']:,.0f}")
print(f"VSS (EEV-RP) = ${VSS:,.0f} ({VSS/tot['EEV']:.1%} of EEV) <- value of the stochastic solution")
print(f"EVPI (RP-WS) = ${EVPI:,.0f} ({EVPI/tot['RP']:.1%}) <- value of perfect information")
fig, ax = plt.subplots(1, 2, figsize=(13, 4.5))
ax[0].bar(["EEV\n(deterministic)","RP\n(stochastic)","WS\n(perfect info)"], tot.values, color=[ORANGE,GREEN,GREY])
for i,v in enumerate(tot.values): ax[0].text(i,v,f"${v:,.0f}",ha="center",va="bottom",fontsize=9)
ax[0].set_ylabel("total realized cost"); ax[0].set_title("Cost by planning approach")
cum = R[["EEV","RP"]].cumsum()
ax[1].plot(range(len(cum)), cum["EEV"], color=ORANGE, lw=2, label="EEV (deterministic)")
ax[1].plot(range(len(cum)), cum["RP"], color=GREEN, lw=2, label="RP (stochastic)")
ax[1].set_xlabel("held-out week #"); ax[1].set_ylabel("cumulative realized cost"); ax[1].set_title("Stochastic plan accumulates less cost"); ax[1].legend(fontsize=8)
fig.tight_layout(); plt.show()
held-out weeks: 41 | RP <= EEV in 39/41 weeks total realized cost EEV $46,521 | RP $44,515 | WS $44,191 VSS (EEV-RP) = $2,006 (4.3% of EEV) <- value of the stochastic solution EVPI (RP-WS) = $324 (0.7%) <- value of perfect information
9 · Does the scenario source matter? — Bayesian vs bootstrap¶
A fair skeptic asks whether the gain comes from stochastic optimization in general, or specifically from the Bayesian scenarios. We repeat the backtest with a frequentist scenario set: bootstrap the historical forecast residuals around the point forecast, instead of drawing from the posterior predictive.
base_tr = np.exp(post["a"].mean()+post["b_dow"].mean(0)[Ftr.dow.values]+post["b_mon"].mean(0)[Ftr.month.values]
+post["b_t"].mean()*Ftr.t.values+post["b_s"].mean(0)[0]*Ftr.sin.values+post["b_s"].mean(0)[1]*Ftr.cos.values)
resid = tr.values - base_tr; rng = np.random.default_rng(1)
def boot_scen(wF, n=200):
base = bayes_mu(wF); return np.clip(base[None,:] + rng.choice(resid, size=(n,7)), 0, None)
tot_boot = np.sum([realized_cost(solve_sp(boot_scen(te_feat.loc[wk.index]), BUDGET), wk.values) for wk in weeks_idx])
vals = [tot['EEV'], tot_boot, tot['RP'], tot['WS']]
print(f"EEV ${tot['EEV']:,.0f} | Bootstrap-RP ${tot_boot:,.0f} | Bayesian-RP ${tot['RP']:,.0f} | WS ${tot['WS']:,.0f}")
fig, ax = plt.subplots(figsize=(8, 4.2))
ax.bar(["EEV\n(mean plan)","Bootstrap RP","Bayesian RP","WS\n(perfect info)"], vals, color=[ORANGE,BLUE,GREEN,GREY])
for i,v in enumerate(vals): ax.text(i,v,f"${v:,.0f}",ha="center",va="bottom",fontsize=8)
ax.set_ylabel("total realized cost"); ax.set_title("Scenario source: both stochastic sets beat the mean plan"); fig.tight_layout(); plt.show()
EEV $46,521 | Bootstrap-RP $44,506 | Bayesian-RP $44,515 | WS $44,191
Honest reading. On this clean, data-rich daily series the Bayesian and bootstrap scenario sets land essentially level — both capture the spread, and with three years of history the posterior is tight enough that it and the empirical residuals largely agree. What clearly matters is being stochastic at all: either scenario set beats the deterministic mean plan by the VSS. The Bayesian route earns its keep where a bootstrap struggles — little history, cold-start units, or extrapolation (a new department, a demand regime shift), where priors and full parameter uncertainty do real work — the same lesson as the predict-then-optimize arc's cold-start finding.
10 · Takeaways¶
- Forecasting is necessary but not sufficient. ARIMA, Holt-Winters, TBATS, a neural net and a Bayesian model all forecast this seasonal series about equally well — a better point forecast was not where the value lay.
- The value is in the decision, and in carrying uncertainty into it. A two-stage stochastic program provisions capacity now and recovers with overtime later; feeding it a distribution rather than a point yields a measurable VSS and beats the naive mean plan in almost every held-out week, while the EVPI bounds what perfect foresight could add.
- Bayesian scenarios integrate naturally into the optimization and match a bootstrap here, with a clear edge reserved for the data-scarce / extrapolation regime.
This is the second family of the ML in Operations Research section: where Predict-then-Optimize put a forecast distribution into a one-shot order, Stochastic Allocation puts a Bayesian scenario set into a two-stage resource decision — the same thesis, one rung up the OR ladder.