C03 · Bike-share demand — Bayesian Modelling in Practice

The brief

The operations team of a bike-sharing scheme plans two things a day ahead: how many vans and staff to put on rebalancing (moving bikes from full docks to empty ones), and how many bikes to have in service. Both scale with the number of rentals they must be able to serve tomorrow. Plan too low and customers find empty docks; plan too high and vans and people stand idle.

They have a weather forecast for tomorrow and two years of history. What they ask for:

"Give us tomorrow's demand with a range we can trust. If you say 90%, we want to be inside it nine days out of ten. And tell us how much capacity to plan so that we can cover demand on 95% of days."

From conversations with the team you know that a typical day sees a few thousand rentals, that a busy summer day can reach eight thousand, and that the system is growing fast.

A point forecast is the easy half of this. The brief is about the spread: a predictive distribution whose stated probabilities can be taken at face value. Sampler diagnostics will not tell you whether you have one.

How this notebook works

Each task states what to deliver, not how. Write your code in the YOUR CODE HERE cells.
Stuck? h.hint("task2") reveals hints one level at a time: nudge → approach → code skeleton. Try to get by on nudges.
h.check("task2", coverage90=...) compares your numbers with the reference solution.
A full worked solution lives in notebooks/solutions/. Open it only when you are done (or truly stuck).

import arviz as az
import matplotlib.pyplot as plt
import numpy as np
import pandas as pd
import pymc as pm

from pymc_challenges import Hints, data

RANDOM_SEED = 2012
rng = np.random.default_rng(RANDOM_SEED)
az.style.use("arviz-variat")

h = Hints("C03")
h.tasks()

C03 - Bike sharing: demand forecasts that admit what they do not know

task0 - Decode, split, look (2 hints, 2 check(s))
task1 - The textbook model for counts (3 hints, 1 check(s))
task2 - Fix what Task 1 revealed (3 hints, 3 check(s))
task3 - Too hot to cycle (3 hints, 2 check(s))
task4 - Calendar, weather situation, growth - and a scoreboard (3 hints, 2 check(s))
task5 - The days the model cannot believe (3 hints, 2 check(s))
task6 - Open the hold-out (3 hints, 3 check(s))
task7 - How much capacity? (3 hints, 3 check(s))

The data

One row per day, 1 January 2011 to 31 December 2012. cnt is the total number of rentals (casual + registered). Calendar columns: yr (0 = 2011, 1 = 2012), mnth, weekday (an integer - find out which day is 0), holiday, workingday (neither weekend nor holiday). weathersit is the day's weather situation: 1 = clear or partly cloudy, 2 = mist or overcast, 3 = light rain or snow, thunderstorm.

The UCI file stores the continuous weather variables normalised - each was divided by a constant so that it lies between 0 and 1:

column	meaning	divided by
`temp`	mean temperature, °C	41
`atemp`	"feels-like" temperature, °C	50
`hum`	relative humidity, %	100
`windspeed`	wind speed, km/h	67

Nobody on the operations team thinks in units of "0.34 of 41 °C".

data.describe("bike_day")
raw = data.load("bike_day")
raw.head()

bike_day: Daily rental counts 2011-2012 with weather and calendar covariates.
  source: Fanaee-T & Gama (2013), Capital Bikeshare (Washington DC), UCI ML Repository.
  url:    https://archive.ics.uci.edu/static/public/275/bike+sharing+dataset.zip

	instant	dteday	season	mnth	weekday	workingday	weathersit	temp	atemp	hum	windspeed	casual	registered	cnt
0	1	2011-01-01	1	1	6	0	2	0.344167	0.363625	0.805833	0.160446	331	654	985
1	2	2011-01-02	1	1	0	0	2	0.363478	0.353739	0.696087	0.248539	131	670	801
2	3	2011-01-03	1	1	1	1	1	0.196364	0.189405	0.437273	0.248309	120	1229	1349
3	4	2011-01-04	1	1	2	1	1	0.200000	0.212122	0.590435	0.160296	108	1454	1562
4	5	2011-01-05	1	1	3	1	1	0.226957	0.229270	0.436957	0.186900	82	1518	1600

The final eight weeks (56 days, from 6 November 2012) are the hold-out. They stay locked away until Task 6: no plots, no fitting, no peeking.

HOLDOUT_START = pd.Timestamp("2012-11-06")

Task 0 · Decode, split, look

Deliver

A data frame with the weather in physical units (°C, %, km/h) and readable labels for weekday and weather situation. Verify the weekday coding rather than guessing it.
train (before HOLDOUT_START) and test (the rest).
Plots of the training data that show what a demand model will have to reproduce: the series over time, demand against temperature, demand by weather situation and by weekday. Write down the three or four features you consider essential. Do you need both temp and atemp?

Your turn Work this task in the notebook.

Hints for this task — open one at a time

1 Nudge

Normalising was a division, so decoding is a multiplication - the constants are in the table above the data. For the weekday coding, do not guess: pandas can tell you the name of the day for every date (.dt.day_name()); cross-tabulate that against weekday.

2 Approach

df.assign(temp_c=df.temp * 41, ...) keeps everything in one tidy frame; wrap it in a function, because in Task 7 you will need to build the same columns for days that are not in the file. pd.Categorical.from_codes(codes, categories) turns integer codes into labels.

For the plots: colour the time series by weather situation, colour the temperature scatter by year, and use box plots for the categorical variables. Look for (1) the shape of the temperature relationship at the hot end, (2) the vertical gap between the years, (3) how much scatter is left at a given temperature, (4) isolated days far below the rest. temp and atemp have a correlation of 0.99: pick one.

Self-checks in the notebook (h.check("task0", …)):

mean_temp_c — mean temperature over all 731 days in degrees Celsius, about 20.3 (temp * 41)
max_wind_kmh — highest daily wind speed in km/h, about 34 (windspeed * 67)

Task 1 · The textbook model for counts

Rentals are counts, so start where every textbook starts: a Poisson regression with a log link. Use the weather only - temperature, humidity and wind speed, each entering linearly on the log scale.

Deliver

Priors you can defend, shown with a prior predictive check on the scale the operations team thinks in: rentals per day.
A fit with clean sampler diagnostics. Look at how precisely the model claims to know its coefficients.
A verdict: is this model's predictive distribution fit for the brief? Back it with at least one graphical posterior predictive check that could embarrass the model, and one number: the share of training days that fall inside the model's central 90% predictive interval.

Your turn Work this task in the notebook.

Hints for this task — open one at a time

1 Nudge

Two separate questions. (1) Did the sampler work? That is what r_hat, ESS and divergences answer. (2) Does the model describe the data? Only a predictive check answers that. A Poisson distribution has a single parameter. Once the mean is fixed, what else is fixed with it - and how wide is a Poisson distribution with a mean of 4,500?

2 Approach

Standardise the three predictors (with the training mean and sd) so that one prior fits all coefficients: on the log scale, a coefficient of 0.5 means "one standard deviation of temperature multiplies demand by e^0.5 = 1.65". Put the intercept prior where the brief puts typical demand: Normal(log(4000), 1). Then pm.sample_prior_predictive and look at quantiles of the simulated counts - try Normal(0, 10) on the intercept as well to see what "vague" means after exp.

For the check, take the posterior predictive draws of cnt, compute the 5% and 95% quantiles per day over chains and draws, and count how many observed days fall in between. Plot that band against the date. az.plot_ppc_pit is the graphical version: if the predictive distributions were right, the PIT values would be uniform and the delta-ECDF curve would hug zero (ArviZ prints the p-value of a uniformity test on the plot). az.plot_ppc_rootogram is made for small counts; with counts in the thousands it is slow and unreadable.

3 Skeleton

Z = (train[cols] - train[cols].mean()) / train[cols].std()
with pm.Model(coords={"predictor": cols}) as poisson_model:
    intercept = pm.Normal("intercept", ___, ___)
    beta = pm.Normal("beta", 0, ___, dims="predictor")
    mu = pm.math.exp(intercept + pm.math.dot(Z.values, beta))
    pm.Poisson("cnt", mu, observed=train.cnt.values)
    prior = pm.sample_prior_predictive(1000, random_seed=RANDOM_SEED)
    idata = pm.sample(random_seed=RANDOM_SEED)
    pm.sample_posterior_predictive(idata, extend_inferencedata=True, random_seed=RANDOM_SEED)

pp = idata.posterior_predictive["cnt"]
q = pp.quantile([0.05, 0.95], dim=("chain", "draw"))
inside = (train.cnt.values >= q.isel(quantile=0).values) & (___)
coverage90 = inside.mean()

Self-checks in the notebook (h.check("task1", …)):

coverage90 — share of training days inside the Poisson model's central 90% predictive interval - about 0.05, nowhere near 0.9

Task 2 · Fix what Task 1 revealed

Your verdict should have been damning, and "the sampler is fine" is no defence.

Deliver

The diagnosis in one sentence: which assumption of the Task 1 model do the data violate, and by roughly what factor? (A number, please.)
A repaired model that changes only what the diagnosis demands - same predictors, same priors for the coefficients. Explain what its new parameter means in words the operations team would understand.
The same checks as in Task 1, now passed (or not?).
Put the posterior standard deviations of the coefficients from both models side by side. What happened, and what would it have cost someone who believed the first model?

Your turn Work this task in the notebook.

Hints for this task — open one at a time

1 Nudge

The Poisson likelihood assumes variance = mean. Compare the two in the data given the model: for each day, how large is (observed - predicted)**2 relative to predicted? If the assumption held, that ratio would average about 1. Which count distribution has a second parameter that lets the variance exceed the mean?

2 Approach

The Pearson dispersion statistic is sum((y - mu)**2 / mu) / (n - number of parameters); use the posterior mean prediction as mu. It should be about 1 for a Poisson model.

The negative binomial has variance mu + mu**2 / alpha. For counts in the thousands the second term dominates, so 1 / sqrt(alpha) is, to a good approximation, the day-to-day coefficient of variation that is left once the predictors are accounted for - a number the operations team can understand ("demand is unpredictable to within about +-x% even if we know the weather"). You can put the prior on that quantity directly and derive alpha from it. pm.NegativeBinomial("cnt", mu=mu, alpha=alpha, observed=...).

Then look at what happened to the posterior sd of the coefficients. The Poisson model treated every one of 3 million rentals as an independent piece of evidence.

3 Skeleton

mu_hat = poisson_idata.posterior_predictive["cnt"].mean(("chain", "draw")).values
dispersion = np.sum((y - mu_hat) ** 2 / mu_hat) / (len(y) - 4)

with pm.Model(coords=coords) as negbin_model:
    intercept = ...   # as before
    beta = ...        # as before
    cv = pm.HalfNormal("cv", ___)
    alpha = pm.Deterministic("alpha", ___)
    mu = pm.math.exp(intercept + pm.math.dot(Z.values, beta))
    pm.NegativeBinomial("cnt", mu=mu, alpha=alpha, observed=train.cnt.values)

Self-checks in the notebook (h.check("task2", …)):

dispersion — Pearson dispersion statistic sum((y - mu)^2 / mu) / (n - p) of the POISSON model from Task 1 - about 500, where 1 would be fine
coverage90 — share of training days inside the repaired model's 90% predictive interval, close to 0.9
alpha — posterior mean of the negative-binomial alpha (variance = mu + mu^2 / alpha) with the three linear weather predictors, about 7

Task 3 · Too hot to cycle

A straight line on the log scale says every extra degree multiplies demand by the same factor, for ever. Washington summers disagree.

Deliver

Evidence from the Task 2 model that the linear temperature effect is wrong (a residual plot will do).
A model in which temperature enters through a smooth function instead. Everything else stays as in Task 2.
A plot of the estimated temperature effect as a multiplier on expected demand, with a credible band, across the observed temperature range.
Two numbers with uncertainty: the temperature at which demand peaks, and how much lower (in %) expected demand is at 35 °C than at the peak.

Your turn Work this task in the notebook.

Hints for this task — open one at a time

1 Nudge

Plot log(observed / predicted) against temperature and average it in bins. If the linear term were adequate, the binned means would hover around zero everywhere. For the fix you need a function of temperature that is flexible but smooth, and that you can evaluate again for new temperatures when you forecast.

2 Approach

A regression spline turns one column into a handful of smooth basis columns; the model stays linear in the coefficients. With patsy: dmatrix("cr(temp_c, df=6, constraints='center') - 1", train) gives a natural cubic spline basis, centred so that it does not compete with the intercept. To evaluate the same basis (same knots) on new data use build_design_matrices([basis.design_info], new_df). Alternative: pm.gp.HSGP with an ExpQuad covariance on standardised temperature.

The temperature effect is f = B_grid @ w on the log scale, so exp(f) is a multiplier on expected demand. Compute it for every posterior draw on a grid of temperatures; then the location of the maximum and the ratio exp(f(35) - max f) are just functions of each draw, and you get their posterior for free.

Register every input with pm.Data now (and give the likelihood shape=mu.shape): you will need to swap in new days in Task 6.

3 Skeleton

from patsy import dmatrix, build_design_matrices
basis = dmatrix("cr(temp_c, df=6, constraints='center') - 1", train)
B = np.asarray(basis)

with pm.Model(coords={"knot": np.arange(B.shape[1]), "day": train.dteday.values}) as m:
    B_ = pm.Data("B", B, dims=("day", "knot"))
    w = pm.Normal("w_temp", 0, 1, dims="knot")
    eta = intercept + pm.math.dot(B_, w) + ___          # humidity and wind as before
    ...

grid = pd.DataFrame({"temp_c": np.linspace(train.temp_c.min(), train.temp_c.max(), 200)})
B_grid = np.asarray(build_design_matrices([basis.design_info], grid)[0])
w_draws = az.extract(idata, var_names="w_temp").transpose("sample", "knot").values
f = w_draws @ B_grid.T                                   # (draws, grid)
peak = grid.temp_c.values[f.argmax(axis=1)]
drop35 = 100 * (1 - np.exp(f[:, ___] - f.max(axis=1)))

Self-checks in the notebook (h.check("task3", …)):

peak_temp_c — posterior mean of the temperature (deg C) at which expected demand peaks, about 26
drop_at_35_pct — posterior mean of how far (in %) expected demand at 35 C lies below the peak, roughly a third

Task 4 · Calendar, weather situation, growth - and a scoreboard

The operations team knows that demand depends on the day of the week, on public holidays, on the weather situation (rain is worse than humidity alone suggests), and that the system grew between 2011 and 2012. Add all four to the Task 3 model.

Deliver

A model in which the categorical effects have no arbitrary reference level: every weekday and every weather situation is treated symmetrically by the prior, and the intercept keeps its meaning. Clean diagnostics. The file also offers workingday - explain why you should not simply add it as well.
The effects in units the team can use: by what percentage does a wet day change demand relative to a clear one? How much did demand grow from 2011 to 2012? Which weekday is busiest? Do not forget the uncertainty.
A PSIS-LOO comparison of all four models so far (az.compare, az.plot_compare) and your reading of it: which steps mattered, judged by elpd_diff against dse? Read the warnings - you will come back to them in Task 5.

Your turn Work this task in the notebook.

Hints for this task — open one at a time

1 Nudge

With dummy variables, one level disappears into the intercept and every other effect is "relative to Sunday". Is there a parameterisation in which seven weekday effects are deviations from the average day? And think about workingday: can you compute it exactly from weekday and holiday? What does a model do with a predictor that is a deterministic function of others?

2 Approach

pm.ZeroSumNormal("weekday_eff", sigma, dims="weekday") is a vector of effects constrained to sum to zero; index it with the integer weekday code: weekday_eff[weekday_idx]. Same for the three weather situations. A contrast such as wet-vs-clear is exp(eff[wet] - eff[clear]) - 1, computed draw by draw.

For growth, the file offers yr. Centre it (-0.5 / +0.5) so that the intercept remains the average day. exp(growth) - 1 is the relative change.

LOO: pm.compute_log_likelihood(idata) inside each model, then az.loo(idata, pointwise=True) once per model and az.compare({name: loo, ...}) - az.compare accepts the precomputed results. A difference matters when elpd_diff is several times its dse. The Poisson row will look absurd, and will come with Pareto-k warnings; both are information, not bugs.

3 Skeleton

coords = {..., "weekday": ["Sun", ...], "weather": ["clear", "mist", "wet"]}
weekday_eff = pm.ZeroSumNormal("weekday_eff", ___, dims="weekday")
weather_eff = pm.ZeroSumNormal("weather_eff", ___, dims="weather")
b_holiday = pm.Normal("b_holiday", 0, 0.5)
growth = pm.Normal("growth", 0, 1)
eta = eta + weekday_eff[weekday_idx] + weather_eff[___] + b_holiday * holiday + growth * (yr - 0.5)

post = az.extract(idata)
wet_vs_clear = 100 * (np.exp(post["weather_eff"].sel(weather="wet") - ___) - 1)

loos = {name: az.loo(idata, pointwise=True) for name, idata in models.items()}
az.compare(loos)

Self-checks in the notebook (h.check("task4", …)):

wet_vs_clear_pct — posterior mean change in demand on a wet day relative to a clear day, in percent (negative, roughly minus 45)
growth_pct — posterior mean growth in demand from 2011 to 2012 at equal weather and calendar, in percent (roughly +50)

Task 5 · The days the model cannot believe

LOO is not only a scoreboard. It computes, for every single day, how well the model predicts that day when it has not seen it - and how much the posterior leans on it.

Deliver

For your best model so far: the pointwise elpd and the Pareto-$k$ value of every training day, plotted against the date, and a table of the six most surprising days.
For the worst of them: what happened? (You have the date and the internet.) Is it the kind of day the brief asks you to forecast? Look at the other five as well - is any of them a data problem rather than a demand problem?
A decision on how to handle what you found, with a justification an auditor would accept, and a refit. What changed in the dispersion and in the width of the predictive intervals - and why does one day out of 675 have that much leverage?

Your turn Work this task in the notebook.

Hints for this task — open one at a time

1 Nudge

az.loo(..., pointwise=True) returns more than a total: elpd_i (one value per day - very negative means "the model found this day almost impossible") and pareto_k (how much the posterior would move if that day were removed). Put both next to the dates. Then ask of each suspicious day: is this demand the model failed to predict, or was the system (or a sensor) not operating normally?

2 Approach

Attach loo.elpd_i.values and loo.pareto_k.values to the training frame and sort. The worst day had 22 rentals. Look at the date and at what happened on the US east coast that day; the day after is affected as well. A humidity of exactly 0% on a rainy day is a different kind of anomaly.

Options: (a) drop the day(s) and state that forecasts are conditional on the system operating; (b) give the day its own indicator (equivalent for prediction); (c) keep it and let the dispersion parameter absorb it. To see what (c) costs, compare cv / alpha and the average width of the 90% predictive interval with and without. One dispersion parameter serves all days - a single "impossible" day can only be accommodated by making every day noisier.

3 Skeleton

loo = az.loo(idata, pointwise=True)
pointwise = train[["dteday", "cnt", "weather", "hum_pct"]].assign(
    elpd=loo.elpd_i.values, pareto_k=loo.pareto_k.values)
pointwise.nsmallest(6, "elpd")

drop = pd.to_datetime([___])
train_clean = train[~train.dteday.isin(drop)].copy()
train_clean["hum_pct"] = train_clean.hum_pct.where(train_clean.hum_pct > 0).interpolate()
# refit the same model on train_clean, compare alpha and interval widths

Self-checks in the notebook (h.check("task5", …)):

cnt_worst_day — number of rentals on the day with the lowest pointwise elpd (an almost empty system)
k_worst_day — Pareto-k of that day, beyond the 0.7 reliability threshold

Task 6 · Open the hold-out

Time to find out whether "90%" means 90%. Forecast the 56 hold-out days with the model from Task 5, using the weather that actually occurred (a perfect weather forecast - the best case).

Deliver

Forecasts produced without refitting: swap the inputs of the fitted model and sample the posterior predictive.
A plot of the 50% and 90% predictive intervals against the actual counts, plus the empirical coverage of both intervals and the mean error (forecast median minus actual).
If the forecasts are biased: find out why - not by tuning against the hold-out, but by asking what your growth term assumes about November 2012, and checking that assumption on the training data. Fix the model, compare old and new with LOO, and forecast again.
A list of the hold-out days that still fall outside the 90% interval. What do they have in common? Could the training data have warned you? (Task 5 again.) Which direction is the error, and how much does that direction hurt the operations team?

Your turn Work this task in the notebook.

Hints for this task — open one at a time

1 Nudge

A term that says "2012 is x% above 2011" describes the past well. What level does it predict for the end of 2012 - and for 2013? If the system grew steadily, where within each year would a constant-per-year model be too high, and where too low? You can check that on the training residuals without touching the hold-out.

2 Approach

Forecasting: inputs registered with pm.Data, likelihood with shape=mu.shape; then pm.set_data({...}, coords={"day": test.dteday.values}) and pm.sample_posterior_predictive(idata, predictions=True) -> pred.predictions["cnt"]. Build the spline basis for the test days with build_design_matrices and standardise with the training means.

Coverage: per-day quantiles over chain and draw, as in Task 1. Mean error: predictive median minus actual.

Diagnosis: average log(observed / predicted) by month on the training data. A saw-tooth (negative early in each year, positive late) means growth happened within the years. Replace the year step by a linear term in continuous time (years since the start, centred); exp(growth) - 1 is then the annual growth rate. Compare the two with az.compare - both are fitted to the same days, so the comparison is fair.

Then list the hold-out days outside the 90% band and look at their dates.

3 Skeleton

def forecast(model, idata, new):
    with model:
        pm.set_data(features(new), coords={"day": new.dteday.values})
        pred = pm.sample_posterior_predictive(idata, predictions=True, random_seed=RANDOM_SEED)
    return pred.predictions["cnt"]

pp = forecast(model, idata, test)
q = pp.quantile([0.05, 0.25, 0.5, 0.75, 0.95], dim=("chain", "draw")).values
y = test.cnt.values
coverage90 = np.mean((y >= q[0]) & (y <= q[4]))
mean_error = np.mean(q[2] - y)

# trend in continuous time instead of a year step
t = (df.dteday - pd.Timestamp("2011-01-01")).dt.days / 365.25 - 1.0
eta = ___ + growth * t

Self-checks in the notebook (h.check("task6", …)):

coverage90_first — hold-out coverage of the 90% interval for the model with a step between years - far too low (about 0.65)
coverage90_final — hold-out coverage of the 90% interval after fixing the growth term (about 0.8; the misses are holiday periods)
growth_pct_per_year — annual growth rate in percent implied by a linear trend in log demand (roughly 60)

Task 7 · How much capacity?

The team plans three days in the coming weeks and hands you the weather forecast for each. (None is a public holiday.)

scenarios = pd.DataFrame(
    {
        "date": pd.to_datetime(["2012-11-16", "2012-11-18", "2012-11-20"]),
        "temp_c": [14.0, 8.0, 11.0],
        "hum_pct": [50.0, 70.0, 90.0],
        "wind_kmh": [10.0, 20.0, 15.0],
        "weather": ["clear", "mist", "wet"],
    },
    index=["mild clear Friday", "cold misty Sunday", "wet Tuesday"],
)
COST_LOST_RENTAL = 4.0  # $ per rental the system could not serve (lost fare and goodwill)
COST_IDLE = 0.5  # $ per unit of capacity planned but not used (vans, staff, bikes standing by)
scenarios

	date	temp_c	hum_pct	wind_kmh	weather
mild clear Friday	2012-11-16	14.0	50.0	10.0	clear
cold misty Sunday	2012-11-18	8.0	70.0	20.0	mist
wet Tuesday	2012-11-20	11.0	90.0	15.0	wet

Deliver

For each scenario the capacity (in rentals per day) that covers demand with probability 95%, from your final model - and, next to it, the number you would get if you (wrongly) used only the uncertainty about the expected demand. How big is the gap?
Management's 95% is a rule of thumb. With the two costs above, find the capacity that minimises expected cost for each scenario, numerically from the predictive draws. Which quantile of the predictive distribution does it correspond to? What ratio of the two costs would make the 95% rule optimal?
Three sentences for the operations team: what to plan, how far to trust it, and when not to.

Your turn Work this task in the notebook.

Hints for this task — open one at a time

1 Nudge

"Cover demand with probability 95%" is a statement about tomorrow's count, not about the model's uncertainty in the expected count. Which distribution describes tomorrow's count? For part 2: if you are one unit short it costs you 4; if you have one unit too many it costs 0.5. At the optimum the expected cost of adding one more unit equals the expected saving.

2 Approach

Build a small data frame for the three scenario days with exactly the columns your feature function needs (weekday code from the date - careful, pandas counts from Monday, the file from Sunday; weather code from the label; time index from the date), then use the Task 6 forecasting function. Capacity for 95% = the 0.95 quantile of the predictive draws for that day. For comparison compute mu from the posterior draws of the parameters and take its 0.95 quantile: that ignores the negative-binomial scatter.

Expected cost of capacity c: mean(4 * max(d - c, 0) + 0.5 * max(c - d, 0)) over predictive draws d. Evaluate on a grid of c and take the minimum. This is the newsvendor problem: the optimum is the quantile cost_short / (cost_short + cost_idle). Invert that formula to find the cost ratio that justifies a 95% rule.

3 Skeleton

new = scenarios.assign(
    dteday=scenarios.date,
    weekday=(scenarios.date.dt.dayofweek + 1) % 7,
    weathersit=[WEATHER.index(w) + 1 for w in scenarios.weather],
    holiday=0, yr=1,
    t=(scenarios.date - pd.Timestamp("2011-01-01")).dt.days / 365.25,
)
pp = forecast(final_model, final_idata, new)
d = pp.stack(sample=("chain", "draw")).transpose("sample", "day").values   # (draws, 3)
capacity95 = np.quantile(d, ___, axis=0)

grid = np.arange(0, 12001, 10)
short = np.clip(d[:, 0][None, :] - grid[:, None], 0, None)
idle = np.clip(___, 0, None)
cost = (COST_LOST_RENTAL * short + COST_IDLE * idle).mean(axis=1)
best = grid[cost.argmin()]

Self-checks in the notebook (h.check("task7", …)):

capacity95_friday — 95% predictive quantile of rentals for the mild clear Friday, about 7,900 (must include the negative-binomial scatter, not just the uncertainty in mu)
optimal_capacity_friday — cost-minimising capacity for the mild clear Friday with costs 4 and 0.5, about 7,300
optimal_quantile — quantile of the predictive distribution at the cost-minimising capacity, 4 / 4.5 = 0.89

Going further

Holiday seasons. The final model over-forecasts Thanksgiving weekend and the week after Christmas. Give it a way to learn that from 2011 (one Christmas is not much - what prior keeps the estimate honest?) and check whether hold-out coverage improves.
Two populations. casual and registered riders have opposite weekday patterns and different sensitivity to weather. Model them separately (or jointly, with shared growth) and add the predictive draws. Does the total forecast get sharper?
A Gaussian process instead of a spline. Replace the temperature spline with pm.gp.HSGP, and the linear trend with a slowly-varying HSGP in time. Does the temperature curve change where data are sparse? What happens to forecasts when the trend is a GP that reverts to its mean?
The weather forecast is uncertain too. Tomorrow's temperature is known to within a couple of degrees and "wet" is a probability. Propagate that through the capacity decision by sampling scenarios. When does it matter?

h.progress()

C03: 0/23 hints used, 18/18 checks passed.

Reference solution

A complete, executed solution with commentary on every task — including the ways the obvious model fails. Try the hint ladder first.

Open the worked solution →

#The brief

#How this notebook works

#The data

#Task 0 · Decode, split, look

#Task 1 · The textbook model for counts

#Task 2 · Fix what Task 1 revealed

#Task 3 · Too hot to cycle

#Task 4 · Calendar, weather situation, growth - and a scoreboard

#Task 5 · The days the model cannot believe

#Task 6 · Open the hold-out

#Task 7 · How much capacity?

#Going further

Reference solution

The brief

How this notebook works

The data

Task 0 · Decode, split, look

Task 1 · The textbook model for counts

Task 2 · Fix what Task 1 revealed

Task 3 · Too hot to cycle

Task 4 · Calendar, weather situation, growth - and a scoreboard

Task 5 · The days the model cannot believe

Task 6 · Open the hold-out

Task 7 · How much capacity?

Going further