Forecast bike-share demand

Case study 1
120 min
Optional · 1 case study needed for the certificate

Case study 1 of 7 · Machine Learning: From Problem to Reliable ModelForecast bike-share demand

Course overview

1. Understand machine learning

2. Frame the prediction problem

3. Evaluate supervised models

4. Build the first models

5. Try other model families

6. Tune and finish the model

7. Understand and choose the model

8. Prepare the model for handoff

9. Case studies

What you will be able to do

Take one real prediction problem from the business question to a tested forecast and a recommendation. You forecast how many bikes a city bike-share system rents in each hour. You also see why the honest answer is less pretty than the first one.

Case studyPortfolio piece

Dataset: Bike Sharing, Washington DC 2011–2012 (UCI)
Licence: CC BY 4.0
Time: About 2 hours
Runs in: Colab, Kaggle or VS Code
Skills used: Leak detection
Time split
Poisson models
Error slices

Open in Colab Open in Kaggle View on GitHub

Colab and Kaggle open your own copy. To use VS Code, download the notebook from GitHub. Stuck? See the solution notebook.

The notebook has every code cell of this page, three tasks with checks, and the instructions. Its setup cells install the packages and download the data, so it needs no local setup. It runs in about one minute. Read the story here, and run the code there.

Before you startWhat to know and to have readyAbout 120 minutes

What you need to know

This case study uses the whole course. The lessons it leans on most:

Why we judge a forecast on later data and use the test set once, from Train, validation and test sets: How to split: at random or by time.
The mean error, MAE and RMSE, and why MAPE breaks near zero, from Regression metrics.
A baseline that a model must beat, from Define a baseline.
Pipelines with one-hot encoding and scaling, from Preprocess with a pipeline.
Gradient boosting, from Ensembles: Boosting: learn from the mistakes, and the experiment log, from Machine learning is iterative.
Leaks, from Select useful features, and slices with counts, from Slice analysis.

What you need on your computer (only for the local route; the notebook needs none of it)

A Python environment with pandas 3.0.6, scikit-learn 1.9.1 and matplotlib, such as the .venv of your project from Define a baseline: Set up your project. Download the notebook from the box above into any folder, and open it as in Define a baseline: Open a notebook in VS Code. The notebook downloads the data itself into a new folder bike-demand/. To start again, delete that folder.

Tested on macOS with Python 3.14.6, pandas 3.0.6, scikit-learn 1.9.1 and matplotlib 3.11.2. The notebook was tested in Jupyter. Windows is not tested.

The situation

Amara Okafor runs operations for a city bike-share system: a few hundred stations where people take a bike and leave it at another station. Every evening she plans the next day. She decides how many staff work at the busiest stations, and how many vans drive around to move bikes from full stations to empty ones.

Both kinds of mistake cost her. If she plans too little for a busy hour, commuters find empty stations at 8 in the morning. Some of them stop using the system. If she plans too much, she pays drivers and vans that wait. The first mistake is worse: a lost customer costs more than an hour of a van. Today she plans from experience and last week’s numbers.

She asks Daniel Kim, the data analyst, a precise question: “How many bikes will people rent in each hour tomorrow?” She needs the answer at 8 in the evening, with tomorrow’s weather forecast. She can afford about two more vans at peak times, so she wants to know when the peaks will be and how big.

What you will build

You will build an hourly demand forecast for the next day, and a rule that turns it into a plan. You will compare it honestly with what a planner can do without machine learning, and you will write down where it fails. At the end you have the parts of a portfolio piece: a notebook, an experiment log with nine experiments, a saved model, a short model card and a recommendation.

The case study uses skills from the whole course. You frame the problem and find a leak (Modules 2 and 6). You split by time and choose regression metrics (Module 3). You build pipelines (Module 4), boosting and an experiment log (Module 5), and tune one setting (Module 6). Then you check slices (Module 7) and write a model card (Module 8). Two ideas are new, and the page explains them when they come: the Poisson loss for counts, and a plan factor that turns a forecast into a decision.

The data

The data is the Bike Sharing dataset from the UCI Machine Learning Repository. It has every hour of the Capital Bikeshare system in Washington, D.C., in 2011 and 2012, with the weather of that hour. Hadi Fanaee-T and João Gama made it from the public trip records of the system. They added the weather from freemeteo.com and the city’s list of holidays.

Credit: Fanaee-T, H. (2013). Bike Sharing [Dataset]. UCI Machine Learning Repository. doi:10.24432/C5W894. Licensed under CC BY 4.0, which allows reuse with credit. The authors ask you to cite their paper: Fanaee-T, H. and Gama, J. (2013), Event labeling combining ensemble detectors and background knowledge, Progress in Artificial Intelligence, doi:10.1007/s13748-013-0040-3. Nextia hosts the file hour.csv unchanged, compressed, so that the notebook gets the same bytes every time.

Amara’s system is not Capital Bikeshare. She and Daniel are invented, to give the data a decision. The numbers are real.

Load it and look at three rows

Every code block on this page is a cell of the notebook. The setup cells have already loaded pandas as pd and numpy as np, and downloaded the file into the working folder.

bikes = pd.read_csv("bike_sharing_hour.csv.gz", parse_dates=["dteday"])
print(bikes.shape)
print(bikes[["dteday", "hr", "workingday", "weathersit", "temp", "casual", "registered", "cnt"]].head(3).to_string())

The output:

(17379, 17)
      dteday  hr  workingday  weathersit  temp  casual  registered  cnt
0 2011-01-01   0           0           1  0.24       3          13   16
1 2011-01-01   1           0           1  0.22       8          32   40
2 2011-01-01   2           0           1  0.22       5          27   32

What this tells us. The file has 17,379 rows and 17 columns. One row is one hour of the whole system: midnight on 1 January 2011 had 16 rentals. The target is cnt, the number of rentals in the hour. casual and registered split it into users without and with a membership. Look at temp: 0.24 is not a temperature in degrees. The authors divided every temperature by 41, so 0.24 means about 10 °C.

The columns we use

Before any model, write down what each column means. These are the columns this case study uses, in plain words:

Column	Meaning
`dteday`, `hr`	the date and the hour (0 to 23)
`weekday`, `workingday`, `holiday`	day of the week (0 = Sunday); 1 if the day is not a weekend or holiday; 1 on a public holiday
`mnth`, `yr`	month (1 to 12); year (0 = 2011, 1 = 2012)
`weathersit`	1 clear or partly cloudy; 2 mist or cloudy; 3 light rain or light snow; 4 heavy rain, ice or snow
`temp`, `atemp`	temperature and “feels like” temperature, scaled (× 41 and × 50 give °C)
`hum`, `windspeed`	humidity (× 100 gives %) and wind speed, scaled
`cnt`	the target: rentals in the hour
`casual`, `registered`	rentals by casual users and by members

We leave out two columns. instant is only the row number, and season repeats what mnth already says.

First checks

Next, the checks that every new dataset needs: empty cells, duplicate rows, gaps, and codes that look wrong.

print("missing values:", int(bikes.isna().sum().sum()))
print("duplicate hours:", int(bikes.duplicated(["dteday", "hr"]).sum()))
calendar_hours = bikes["dteday"].nunique() * 24
print("hours in the calendar:", calendar_hours, "| rows:", len(bikes), "| missing hours:", calendar_hours - len(bikes))
print("weather codes:", bikes["weathersit"].value_counts().sort_index().to_dict())
print("humidity 0:", int((bikes["hum"] == 0).sum()), "| wind speed 0:", int((bikes["windspeed"] == 0).sum()))

The output:

missing values: 0
duplicate hours: 0
hours in the calendar: 17544 | rows: 17379 | missing hours: 165
weather codes: {1: 11413, 2: 4544, 3: 1419, 4: 3}
humidity 0: 22 | wind speed 0: 2180

What this tells us. No cell is empty and no hour appears twice. But two years have 17,544 hours, and the file has 17,379 rows, so 165 hours have no row at all. Weather code 4 appears in only 3 hours, too few for any model to learn from. A humidity of 0 is not possible in Washington, so those 22 values are sensor gaps. A wind speed of exactly 0 in 2,180 hours is suspicious too: it probably means “not measured”. We keep these rows, and we remember that wind speed is not reliable.

Where are the missing hours? Count them by day:

times = bikes["dteday"] + pd.to_timedelta(bikes["hr"], unit="h")
all_hours = pd.date_range(times.min(), times.max(), freq="h")
missing = pd.Series(all_hours.difference(times))
print(missing.dt.date.value_counts().head(5).to_string())

The output:

2012-10-29    23
2011-01-27    16
2012-10-30    13
2011-01-18    12
2011-01-26     8

What this tells us. The largest gaps are on 29 and 30 October 2012. That was Hurricane Sandy: the system closed. The other gaps, in January 2011, are snow days and data problems. So the data has no row for an hour when nobody could rent a bike. Those two Sandy days are in the last quarter of 2012, which will be our test period.

The target, and a leak

Now look at the target, and check one suspicion: do casual and registered add up to cnt?

print(bikes["cnt"].describe().round(1).to_dict())
print("casual + registered == cnt in every row:", bool((bikes["casual"] + bikes["registered"] == bikes["cnt"]).all()))

The output:

{'count': 17379.0, 'mean': 189.5, 'std': 181.4, 'min': 1.0, '25%': 40.0, '50%': 142.0, '75%': 281.0, 'max': 977.0}
casual + registered == cnt in every row: True

What this tells us. A typical hour has about 142 rentals, and the busiest hour has 977. The smallest count is 1, never 0: an hour with no rentals probably has no row at all. And casual + registered equals the target in every row. These two columns are counted in the same hour as the target, so at 8 in the evening, they are not known for tomorrow. They are a leak: information that is not available at the prediction time, and that makes a model look much better than it can be. You met leaks in Select useful features. They never go into the inputs.

Frame the problem

Now write the problem down in the terms of Translate a need, before any model:

Decision: how many staff and vans Amara plans for each hour of tomorrow.
Prediction unit: one hour of the whole system.
Target: cnt, the rentals in that hour. It is a number, so this is regression.
Prediction time: 20:00 the evening before. Everything an input uses must be known then.
Inputs available at that time: the calendar (hour, weekday, month, year, holiday, working day) and the weather forecast for tomorrow.
Not available: casual and registered (the leak), and the rentals of tomorrow’s earlier hours.
People: no column describes a person, so there is no fairness slice by group. The slices that matter are hours, weather and special days.

One honest limit appears already. The data has the weather that really happened, not the forecast that Amara will have at 20:00. A weather forecast for tomorrow is usually close, but it can miss a storm. So every error on this page is a little better than the error Amara will see. We write that down now, and again in the limits.

Explore what the decision needs

We do not explore every column. We ask three questions that Amara’s plan depends on. What does a day look like? Is the demand changing? How much does the weather matter?

The shape of a day

Amara plans hour by hour, so the first question is the shape of a day. Compute the mean rentals for each hour, separately for working days and for other days:

bikes["day_type"] = np.where(bikes["workingday"] == 1, "working day", "other day")
by_hour = bikes.pivot_table(index="hr", columns="day_type", values="cnt", aggfunc="mean").round(0)
print(by_hour.loc[[4, 8, 13, 17, 22]].to_string())

The output (five of the 24 hours):

day_type  other day  working day
hr
4               8.0          5.0
8             106.0        477.0
13            373.0        198.0
17            324.0        525.0
22            116.0        138.0

The chart shows all 24 hours.

Working day
Weekend or holiday

Mean rentals for each hour of the day, 2011 and 2012. Working days have two commuter peaks; other days have one wide hill around midday.

Text description of the chart

A line chart. The horizontal axis is the hour of the day, from 0 to 23. The vertical axis is the mean rentals per hour, from 0 to 600. The working-day line is low at night, about 5 rentals at 3:00 and 4:00. It rises sharply to a peak of 477 at 8:00 and falls to about 135 to 200 in the middle of the day. It rises again to its highest point, 525 at 17:00 and 492 at 18:00, then falls to 89 at 23:00. The line for weekends and holidays is low until 7:00. It climbs to a wide hill of about 355 to 373 between 12:00 and 16:00, and falls to 86 at 23:00. At 8:00, a working day has about four and a half times the rentals of another day. At 13:00, other days have almost twice as many.

What this tells us. A working day has two peaks, at 8:00 and at 17:00 to 18:00: people ride to work and back. Weekends and holidays have one wide hill around midday: people ride for fun. So the hour alone is not enough. A forecast must know which kind of day it is. It must also get the two commuter peaks right, because those are the hours when stations run empty.

The trend

Next, is the demand changing over time? Compare the mean rentals per hour in the two years:

print("mean rentals per hour:", bikes.groupby(bikes["dteday"].dt.year)["cnt"].mean().round(1).to_dict())
monthly = bikes.groupby(bikes["dteday"].dt.to_period("M"))["cnt"].sum()
print("quietest month:", monthly.idxmin(), int(monthly.min()), "| busiest month:", monthly.idxmax(), int(monthly.max()))

The output:

mean rentals per hour: {2011: 143.8, 2012: 234.7}
quietest month: 2011-01 38189 | busiest month: 2012-09 218573

The mean grew from 143.8 to 234.7 rentals per hour, by 63% in one year. The chart shows the mean rentals per day for each week, coloured by the part of the split that the next section makes.

Train
Valid
Test

Mean rentals per day, week by week, January 2011 to December 2012. Each summer is busier than each winter, and 2012 is busier than 2011 in every season.

Text description of the chart

A line chart of 106 weeks, from the end of December 2010 to the end of December 2012. The vertical axis is the mean rentals per day, from 0 to 8,000. The train part, to June 2012, starts below 1,000 rentals a day in January 2011. It rises to about 4,000 to 5,000 in the summer of 2011 and falls to about 3,000 in the winter. Then it climbs again to about 7,000 in June 2012. The valid part, July to September 2012, is the highest of all: between 6,036 and 7,860 rentals a day. The test part, October to December 2012, starts at about 6,500 to 7,000 and drops to 4,109 in the week of Hurricane Sandy. It stays between about 3,900 and 5,800 in November and early December, and falls to 1,531 in the week of Christmas.

What this tells us. Two patterns overlap: the seasons, and strong growth. Summer 2012 is much busier than summer 2011, because the system added stations and members. This matters for every model. A model learns the level of the past, so it will tend to be too low in the future. And the last quarter of 2012 is both a busy year and a falling season.

The weather

Last, how much does the weather change the demand?

weather = bikes.groupby("weathersit")["cnt"].agg(hours="size", mean_rentals="mean").round(1)
print(weather.to_string())

The output:

            hours  mean_rentals
weathersit
1           11413         204.9
2            4544         175.2
3            1419         111.6
4               3          74.3

What this tells us. In light rain or snow (code 3), the system rents about 112 bikes per hour, almost half of the 205 in clear weather. So the weather forecast for tomorrow is an important input. Rainy hours may be hard to forecast, because rain often starts and stops within an hour.

Split and baseline

A split by time

In Train, validation and test sets, you split by time when the model will predict the future. That is exactly Amara’s case: on any evening, the model knows the past and forecasts tomorrow. So the split follows the calendar:

Train: from January 2011 to the end of June 2012. The models learn from it.
Valid: July to September 2012. We use it to compare models and make every decision.
Test: October to December 2012. It stays sealed until the final test.

TARGET = "cnt"
FEATURES = ["hr", "weekday", "workingday", "holiday", "mnth", "yr", "weathersit", "temp", "atemp", "hum", "windspeed"]
train = bikes[bikes["dteday"] < "2012-07-01"]
valid = bikes[(bikes["dteday"] >= "2012-07-01") & (bikes["dteday"] < "2012-10-01")]
test = bikes[bikes["dteday"] >= "2012-10-01"]
print("train", len(train), "| valid", len(valid), "| test", len(test))
print("mean per hour: train", round(train[TARGET].mean(), 1), "| valid", round(valid[TARGET].mean(), 1))

The output:

train 13003 | valid 2208 | test 2168
mean per hour: train 167.7 | valid 288.4

What this tells us. Train has 13,003 hours (18 months), and valid and test have about 2,200 each (3 months). Valid is much busier than train: 288.4 rentals per hour against 167.7. That is summer plus growth. A model that learned mostly from quieter months will be tested on the busiest quarter of all.

How we measure errors

Which metric fits Amara’s decision? She reads errors in bikes, so the main number is the MAE, the typical miss in bikes per hour. A forecast that is too low leaves stations empty, so the direction matters. We always report the mean error next to it: positive means the forecast is too low. Large misses at the peaks are the worst, so we also report the RMSE, which punishes them. We do not use MAPE: at night, an hour can have 2 rentals, and MAPE breaks near zero, as you saw in Regression metrics.

from sklearn.metrics import mean_absolute_error, root_mean_squared_error

def report(real, forecast):
    """Mean error (positive = too low), MAE and RMSE, in bikes per hour."""
    return {"mean error": round(float((real - forecast).mean()), 1),
            "MAE": round(float(mean_absolute_error(real, forecast)), 1),
            "RMSE": round(float(root_mean_squared_error(real, forecast)), 1)}

Two baselines

A model is worth building only if it beats what Amara can do without it. You learned this in Define a baseline. Daniel tries two baselines. The first is the mean of train, for every hour: the floor. The second is what a careful planner would do in a spreadsheet: the hour profile. It is the mean rentals for each hour of a working day and of other days, over the last 8 weeks of train. Recent weeks follow the growth better than the whole history.

mean_forecast = np.full(len(valid), train[TARGET].mean())
recent = train[train["dteday"] >= "2012-05-06"]  # the last 8 weeks of train
profile = recent.groupby(["workingday", "hr"])[TARGET].mean()
profile_forecast = profile.reindex(pd.MultiIndex.from_frame(valid[["workingday", "hr"]])).to_numpy()
mean_result = report(valid[TARGET], mean_forecast)
profile_result = report(valid[TARGET], profile_forecast)
print("mean of train:", mean_result)
print("hour profile: ", profile_result)

The output:

mean of train: {'mean error': 120.6, 'MAE': 197.8, 'RMSE': 260.7}
hour profile:  {'mean error': 14.8, 'MAE': 52.2, 'RMSE': 79.8}

What this tells us. The mean of train misses by about 198 bikes in a typical hour: it ignores the shape of the day. The hour profile misses by about 52 bikes per hour. It is too low by about 15 on average, because the summer kept growing after the 8 weeks. The profile is a strong baseline, and it costs nothing. A model must beat 52.2, not 197.8.

The experiment log

From here on, every attempt goes into an experiment log, as in Machine learning is iterative. The columns are the same as there, with the metrics of this problem: id, date, change, valid_bias (the mean error), valid_mae, valid_rmse, decision and reason. With data in time order, the validation quarter does the job that cross-validation did for the tickets.

import datetime
from pathlib import Path

LOG_FILE = Path("experiment_log.csv")
LOG_COLUMNS = ["id", "date", "change", "valid_bias", "valid_mae", "valid_rmse", "decision", "reason"]

def log_experiment(id, change, result, decision, reason):
    """Add one row to the log, or replace the row with the same id."""
    log = pd.read_csv(LOG_FILE) if LOG_FILE.exists() else pd.DataFrame(columns=LOG_COLUMNS)
    row = {"id": id, "date": datetime.date.today().isoformat(), "change": change,
           "valid_bias": result["mean error"], "valid_mae": result["MAE"], "valid_rmse": result["RMSE"],
           "decision": decision, "reason": reason}
    log = pd.concat([log[log["id"] != id], pd.DataFrame([row])], ignore_index=True)
    log.sort_values("id").to_csv(LOG_FILE, index=False)

def show_log():
    log = pd.read_csv(LOG_FILE)
    print(log[["id", "change", "valid_bias", "valid_mae", "valid_rmse", "decision"]].to_string(index=False))

LOG_FILE.unlink(missing_ok=True)  # start a new log

log_experiment("E01", "Baseline: the mean of train for every hour", mean_result, "baseline", "The floor.")
log_experiment("E02", "Baseline: mean per hour and day type, last 8 weeks of train", profile_result, "baseline",
               "What a planner can do with a spreadsheet; the model must beat it.")
show_log()

The output:

 id                                                      change  valid_bias  valid_mae  valid_rmse decision
E01                  Baseline: the mean of train for every hour       120.6      197.8       260.7 baseline
E02 Baseline: mean per hour and day type, last 8 weeks of train        14.8       52.2        79.8 baseline

What this tells us. The log now holds the two baselines, with the reason for each. Every model from here on must beat E02, with a valid MAE below 52.2.

Models

Daniel tries three kinds of model, from simple to complex, and keeps every result. Each model is a scikit-learn pipeline: the preprocessing and the model in one object, fitted on train only.

from sklearn.compose import ColumnTransformer
from sklearn.ensemble import HistGradientBoostingRegressor
from sklearn.linear_model import PoissonRegressor, Ridge
from sklearn.pipeline import Pipeline
from sklearn.preprocessing import OneHotEncoder, StandardScaler

CATEGORIES = ["hr", "weekday", "mnth", "weathersit"]
NUMBERS = ["workingday", "holiday", "yr", "temp", "atemp", "hum", "windspeed"]

def linear_pipeline(model, categories=CATEGORIES):
    """One-hot codes for the categories, scaled numbers, then a linear model."""
    prepare = ColumnTransformer([("categories", OneHotEncoder(handle_unknown="ignore"), categories),
                                 ("numbers", StandardScaler(), NUMBERS)])
    return Pipeline([("prepare", prepare), ("model", model)])

def boosting_pipeline(model, features=FEATURES):
    """Keep only the listed columns, then the model. Trees need no scaling."""
    select = ColumnTransformer([("features", "passthrough", features)])
    return Pipeline([("select", select), ("model", model)])

Both functions choose their own columns. linear_pipeline() one-hot encodes the hour, weekday, month and weather code, and scales the numbers, as in Preprocess with a pipeline. boosting_pipeline() only keeps the listed columns, because trees need no scaling. Because each pipeline selects its own inputs, you can pass the whole DataFrame, and casual or registered can never get in by mistake.

A linear model, and why it fails

The first model is ridge regression, a linear regression with a small penalty on large coefficients (you met the penalty as regularization in Control complexity).

pipeline = linear_pipeline(Ridge(alpha=1.0)).fit(train, train[TARGET])
ridge_forecast = pipeline.predict(valid)
ridge_result = report(valid[TARGET], ridge_forecast)
print(ridge_result)

The output:

{'mean error': 17.5, 'MAE': 111.4, 'RMSE': 146.1}

What this tells us. The MAE is 111.4: twice as bad as the hour profile. Why? Compare the forecast with the real mean for four hours of a working day:

check = valid.assign(forecast=ridge_forecast)
hours = check[check["workingday"] == 1].groupby("hr")[[TARGET, "forecast"]].mean().round(0)
print(hours.loc[[3, 8, 13, 17]].to_string())

The output:

      cnt  forecast
hr
3     7.0     109.0
8   693.0     425.0
13  294.0     341.0
17  813.0     518.0

What this tells us. At 3:00, about 7 people rent a bike, and the model says 109. At 17:00, 813 rent a bike, and it says 518. The model is too flat: too high at night, too low at the peaks. The reason is the way a linear model works. It adds effects: “a working day adds so many bikes, rain takes away so many bikes, 8:00 adds so many”. But demand does not work like that.

Counts multiply: the Poisson model

Think of an ice-cream shop. On a hot day, it sells twice as much as usual, in the morning and in the afternoon. It does not sell “200 more” in every hour. At 7:00, when it usually sells 10, it sells 20. At 15:00, when it usually sells 300, it sells 600. Weather and season multiply the demand of each hour. Bikes are the same: rain takes about half of the rentals, in a quiet hour and in a busy one.

A model where effects multiply is a linear model on the logarithm of the expected count. On the log scale, “× 2” becomes “+ log 2”, so adding becomes multiplying again. The model for counts that does this is Poisson regression. It is named after Siméon Poisson, a French mathematician who studied counts of events, such as the number of letters lost in the post each day. Its forecasts are also never below zero, which is right for a count.

pipeline = linear_pipeline(PoissonRegressor(alpha=1e-4, max_iter=1000)).fit(train, train[TARGET])
poisson_result = report(valid[TARGET], pipeline.predict(valid))
print(poisson_result)

The output:

{'mean error': -28.0, 'MAE': 93.0, 'RMSE': 131.9}

What this tells us. The MAE falls from 111.4 to 93.0. Better, but still far from the profile’s 52.2, and now the forecast is too high on average (mean error −28.0). The model still uses one daily shape for every day. It cannot know that 8:00 is a peak on a Monday and quiet on a Sunday. Daniel gives it that information directly, as one category for each hour of each kind of day, 48 values in all:

def add_hour_type(part):
    """One category per hour of a working day and per hour of other days: 48 values."""
    return part.assign(hour_type=part["workingday"].astype(str) + "-" + part["hr"].astype(str))

train, valid, test = add_hour_type(train), add_hour_type(valid), add_hour_type(test)
pipeline = linear_pipeline(PoissonRegressor(alpha=1e-4, max_iter=1000), CATEGORIES + ["hour_type"])
pipeline.fit(train, train[TARGET])
hour_type_result = report(valid[TARGET], pipeline.predict(valid))
print(hour_type_result)

The output:

{'mean error': -27.9, 'MAE': 54.8, 'RMSE': 87.3}

What this tells us. The MAE falls to 54.8. The shape of the day was the missing piece. A combination of two columns that acts differently together is called an interaction. Linear models need you to build it; we will see that trees find it alone. Even so, the best linear model is still a little worse than the hour profile. Into the log:

log_experiment("E03", "Ridge, one-hot hour, weekday, month, weather", ridge_result, "reject",
               "Far worse than the hour profile: it adds effects, but rentals multiply.")
log_experiment("E04", "Poisson regression, the same inputs", poisson_result, "reject",
               "Better than ridge, still far worse than the hour profile.")
log_experiment("E05", "E04 + hour_type (hour x working day)", hour_type_result, "reject",
               "The shape of the day was the missing piece; still worse than the profile.")
show_log()

The output:

 id                                                      change  valid_bias  valid_mae  valid_rmse decision
E01                  Baseline: the mean of train for every hour       120.6      197.8       260.7 baseline
E02 Baseline: mean per hour and day type, last 8 weeks of train        14.8       52.2        79.8 baseline
E03                Ridge, one-hot hour, weekday, month, weather        17.5      111.4       146.1   reject
E04                         Poisson regression, the same inputs       -28.0       93.0       131.9   reject
E05                        E04 + hour_type (hour x working day)       -27.9       54.8        87.3   reject

What this tells us. Three linear models, and none of them beats the hour profile. The log keeps all three, with the reason, so that nobody tries them again without a new idea.

Gradient boosting, with the right loss

Next, gradient boosting, from Ensembles: many small trees, where each new tree corrects the errors of the trees before it. A tree can split first on the kind of day and then on the hour, so it finds the interaction without help. scikit-learn’s HistGradientBoostingRegressor is the regression version of the classifier you used.

Boosting measures its errors with a loss, as linear regression did in Regression intuition. The default is the squared error. Daniel also tries loss="poisson": the same idea as Poisson regression. It judges a miss relative to the size of the count: a miss of 20 bikes in a quiet hour counts more than a miss of 20 in a busy hour. Its forecasts are never below zero.

boosting = boosting_pipeline(HistGradientBoostingRegressor(random_state=0)).fit(train, train[TARGET])
forecast = boosting.predict(valid)
squared_result = report(valid[TARGET], forecast)
print("squared error:", squared_result, "| hours below 0:", int((forecast < 0).sum()))
boosting = boosting_pipeline(HistGradientBoostingRegressor(loss="poisson", random_state=0)).fit(train, train[TARGET])
forecast = boosting.predict(valid)
boosting_result = report(valid[TARGET], forecast)
print("Poisson loss: ", boosting_result, "| hours below 0:", int((forecast < 0).sum()))

The output:

squared error: {'mean error': 20.1, 'MAE': 40.9, 'RMSE': 60.8} | hours below 0: 10
Poisson loss:  {'mean error': 21.4, 'MAE': 40.3, 'RMSE': 62.4} | hours below 0: 0

What this tells us. Both boosting models beat the hour profile by about 12 bikes per hour: 40.9 and 40.3 against 52.2. The squared-error model forecasts fewer than 0 bikes in 10 night hours. That is impossible, and a planner who sees it stops trusting the tool. The Poisson loss never does that, and its MAE is a little lower, so Daniel keeps it. Both models are too low by about 20 bikes per hour on average: valid is busier than anything they learned from.

log_experiment("E06", "Gradient boosting, squared error, defaults", squared_result, "reject",
               "Good, but it forecasts fewer than 0 bikes in some hours.")
log_experiment("E07", "Gradient boosting, Poisson loss, defaults", boosting_result, "keep",
               "Best MAE so far, never below 0; 12 bikes better than the profile.")
show_log()

The output:

 id                                                      change  valid_bias  valid_mae  valid_rmse decision
E01                  Baseline: the mean of train for every hour       120.6      197.8       260.7 baseline
E02 Baseline: mean per hour and day type, last 8 weeks of train        14.8       52.2        79.8 baseline
E03                Ridge, one-hot hour, weekday, month, weather        17.5      111.4       146.1   reject
E04                         Poisson regression, the same inputs       -28.0       93.0       131.9   reject
E05                        E04 + hour_type (hour x working day)       -27.9       54.8        87.3   reject
E06                  Gradient boosting, squared error, defaults        20.1       40.9        60.8   reject
E07                   Gradient boosting, Poisson loss, defaults        21.4       40.3        62.4     keep

What this tells us. E07 is the first row marked keep: the first model that beats the hour profile, by 12 bikes per hour.

Your turn: past rentals as inputs

Daniel has one more idea. The rentals of the recent past should say something about tomorrow. Which past rentals are known at 20:00? The same hour yesterday is not. At 20:00 today, the hours from 20:00 to midnight have not happened yet, and tomorrow’s 21:00 would need today’s 21:00. The same hour 2 days before and 7 days before are known. So he adds both, with a small helper:

def add_past_rentals(part, days):
    """Rentals in the same hour `days` days earlier (empty if that hour is missing)."""
    times = bikes["dteday"] + pd.to_timedelta(bikes["hr"], unit="h")
    rentals = pd.Series(bikes["cnt"].to_numpy(), index=times)
    part_times = part["dteday"] + pd.to_timedelta(part["hr"], unit="h")
    past = rentals.reindex(part_times - pd.Timedelta(days=days)).to_numpy()
    return part.assign(**{f"rentals_{days}d_before": past})

In the notebook, you add the two columns to train and valid, fit the Poisson boosting model with them, and log the result as E08. Predict first: will the past rentals help?

Show the answer

with_past = [add_past_rentals(add_past_rentals(part, 2), 7) for part in (train, valid)]
past_features = FEATURES + ["rentals_2d_before", "rentals_7d_before"]
model = HistGradientBoostingRegressor(loss="poisson", random_state=0)
pipeline = boosting_pipeline(model, past_features).fit(with_past[0], with_past[0][TARGET])
past_result = report(with_past[1][TARGET], pipeline.predict(with_past[1]))
print(past_result)
log_experiment("E08", "E07 + rentals in the same hour 2 and 7 days before", past_result, "reject",
               "Worse than E07: the weather of 2 or 7 days ago adds noise.")

The output:

{'mean error': 19.5, 'MAE': 45.3, 'RMSE': 72.7}

What this tells us. The past rentals make the forecast worse: MAE 45.3 against 40.3. The rentals of two days ago carry that day’s weather, a rainy Tuesday or a sunny Saturday, which tomorrow does not share. The model gets an input that sounds useful but is mostly noise for this decision. A rejected idea is still a result, so it stays in the log. You would not know this without trying, and you would not trust it without the log.

Tune one setting, and check it twice

The last experiment slows down the learning: a learning rate of 0.05 and 400 trees, instead of 0.1 and at most 100. In Tune hyperparameters you compared settings with five folds. Here the data is in time order, and one validation quarter is one noisy result. So Daniel checks the setting on a second quarter too: learn from the data before April 2012, and forecast April to June 2012. This is called a backtest: a test on a past period, as if you had been there.

slower = HistGradientBoostingRegressor(loss="poisson", learning_rate=0.05, max_iter=400, random_state=0)
slower_result = report(valid[TARGET], boosting_pipeline(slower).fit(train, train[TARGET]).predict(valid))
print("E09 on valid:", slower_result)
# Check the choice on an earlier quarter too: learn before April 2012, forecast April-June 2012.
early, spring = train[train["dteday"] < "2012-04-01"], train[train["dteday"] >= "2012-04-01"]
for name, model in [("E07", HistGradientBoostingRegressor(loss="poisson", random_state=0)), ("E09", slower)]:
    print(name, "on April-June 2012:", report(spring[TARGET], boosting_pipeline(model).fit(early, early[TARGET]).predict(spring)))

The output:

E09 on valid: {'mean error': 21.6, 'MAE': 38.7, 'RMSE': 61.0}
E07 on April-June 2012: {'mean error': 13.2, 'MAE': 39.3, 'RMSE': 61.8}
E09 on April-June 2012: {'mean error': 8.6, 'MAE': 37.6, 'RMSE': 59.4}

What this tells us. The slower model wins on both quarters: 38.7 against 40.3 on valid, and 37.6 against 39.3 in the spring. The gains are small, but they agree, so E09 is the final setting.

log_experiment("E09", "E07 with learning_rate 0.05 and 400 trees", slower_result, "keep",
               "Better than E07 on valid and on April-June 2012; the final settings.")
show_log()

The output:

 id                                                      change  valid_bias  valid_mae  valid_rmse decision
E01                  Baseline: the mean of train for every hour       120.6      197.8       260.7 baseline
E02 Baseline: mean per hour and day type, last 8 weeks of train        14.8       52.2        79.8 baseline
E03                Ridge, one-hot hour, weekday, month, weather        17.5      111.4       146.1   reject
E04                         Poisson regression, the same inputs       -28.0       93.0       131.9   reject
E05                        E04 + hour_type (hour x working day)       -27.9       54.8        87.3   reject
E06                  Gradient boosting, squared error, defaults        20.1       40.9        60.8   reject
E07                   Gradient boosting, Poisson loss, defaults        21.4       40.3        62.4     keep
E08          E07 + rentals in the same hour 2 and 7 days before        19.5       45.3        72.7   reject
E09                   E07 with learning_rate 0.05 and 400 trees        21.6       38.7        61.0     keep

The chart summarises the log. E01, the mean of train, is left out: its MAE of 197.8 would squash all the other bars.

Valid MAE of each experiment from E02 to E09. Grey: the hour-profile baseline. Blue: models. Orange: the final setting, E09.

Text description of the chart

A column chart of eight experiments, with the valid MAE in bikes per hour from 0 to 120. E02, the hour profile, the baseline: 52.2. E03, ridge regression: 111.4, the tallest bar. E04, Poisson regression: 93.0. E05, Poisson regression with the hour type: 54.8. E06, boosting with squared error: 40.9. E07, boosting with the Poisson loss: 40.3. E08, E07 with the past rentals: 45.3. E09, the slower boosting: 38.7, the shortest bar. Only the four boosting models are below the baseline.

What this tells us. The log tells the story of the work in one picture. The linear models only came close to the baseline once they knew the shape of each kind of day. Every boosting model beat it, and the two settings marked keep are the two steps that helped.

Common mistake

The random split that looks great

Symptom: Lena, a new analyst on the team, follows a popular tutorial. She shuffles all the hours and keeps 20% for testing. Her model is off by 27 bikes per hour, much better than Daniel’s 38.7. To do even better, she adds casual, and the error falls to 21.

from sklearn.model_selection import train_test_split

shuffled_train, shuffled_test = train_test_split(bikes, test_size=0.2, random_state=0)
model = HistGradientBoostingRegressor(loss="poisson", random_state=0)
pipeline = boosting_pipeline(model).fit(shuffled_train, shuffled_train[TARGET])
print("random split:", report(shuffled_test[TARGET], pipeline.predict(shuffled_test)))
pipeline = boosting_pipeline(model, FEATURES + ["casual"]).fit(shuffled_train, shuffled_train[TARGET])
print("random split + casual:", report(shuffled_test[TARGET], pipeline.predict(shuffled_test)))

The output:

random split: {'mean error': 0.7, 'MAE': 27.3, 'RMSE': 45.6}
random split + casual: {'mean error': 0.8, 'MAE': 21.4, 'RMSE': 34.7}

Cause: with a random split, every test hour has its neighbours in train: the hour before, the hour after, the same hour yesterday. The model learns the level of the very same days, including the growth of 2012 and every storm. Amara never has that: at 20:00, tomorrow’s neighbours do not exist yet. With casual, the model also sees part of the answer. Both numbers measure a situation that will never happen.

Fix: split by time, as the model will be used: learn from the past and test on a later period. Remove every input that is not known at the prediction time. Then expect a larger, honest error: here 38.7 on valid, and more on test.

Consequence in AI work: a forecast that was tested with a random split looks excellent in the report and disappoints in its first month. The team then blames the model, the data or the users, when the real cause was the test.

Metric and operating point

The chosen model, E09, fitted on train and checked on valid:

chosen = HistGradientBoostingRegressor(loss="poisson", learning_rate=0.05, max_iter=400, random_state=0)
pipeline = boosting_pipeline(chosen).fit(train, train[TARGET])
valid_forecast = pipeline.predict(valid)
print(report(valid[TARGET], valid_forecast))

The output:

{'mean error': 21.6, 'MAE': 38.7, 'RMSE': 61.0}

What this tells us. These are the numbers of E09 in the log. The typical miss is 38.7 bikes per hour, and the forecast is too low by 21.6 on average.

A forecast is not yet a plan. The model is too low on average, and for Amara a shortage is worse than a spare bike. So she does not plan with the forecast itself. Think of a baker who knows that the shop sells about 100 loaves on a Saturday. She bakes 120, because running out loses customers, and a few spare loaves cost little. The extra share is her safety margin.

Amara does the same with a plan factor. In the busy hours, the hours with a forecast of 100 bikes or more, she plans for forecast × factor. Daniel chooses the factor on valid, so that the plan covers 9 of 10 busy hours. This is the regression version of the threshold in Thresholds and calibration. It is one number that turns the model’s output into a decision, chosen on valid from the cost of each error. Such a number is the operating point.

busy = valid_forecast >= 100
ratio = valid[TARGET].to_numpy()[busy] / valid_forecast[busy]
plan_factor = round(float(np.quantile(ratio, 0.9)), 2)
print("busy hours in valid:", int(busy.sum()), "| plan factor:", plan_factor)
covered = (valid[TARGET].to_numpy()[busy] <= plan_factor * valid_forecast[busy]).mean()
print("busy hours covered by forecast x factor:", round(float(covered), 3))

The output:

busy hours in valid: 1618 | plan factor: 1.26
busy hours covered by forecast x factor: 0.891

What this tells us. With a factor of 1.26, the plan covers 89.1% of the 1,618 busy hours of valid. In a busy valid hour, that means about 92 bikes of extra capacity on average. That is the price of the safety margin, and Amara can compare it with the vans she has.

Errors and slices

Where does the model fail? As in Slice analysis, we look at groups of hours on valid, with the count of hours next to each error. Every decision uses valid; the test slices come later, and only as a report.

RUSH = [7, 8, 9, 16, 17, 18, 19]

def time_of_day(part):
    rush = (part["workingday"] == 1) & part["hr"].isin(RUSH)
    day = part["hr"].between(7, 19)
    return np.select([rush, day], ["rush hour", "other day hour"], "night (20-6)")

def slice_table(part, forecast, by):
    errors = part.assign(error=part[TARGET] - forecast, abs_error=(part[TARGET] - forecast).abs())
    return errors.groupby(by).agg(hours=("error", "size"), mean_rentals=(TARGET, "mean"),
                                  MAE=("abs_error", "mean"), mean_error=("error", "mean")).round(1)

print(slice_table(valid, valid_forecast, time_of_day(valid)).to_string())

The output:

                hours  mean_rentals   MAE  mean_error
night (20-6)     1012         123.3  20.6        11.0
other day hour    755         335.8  41.8        14.8
rush hour         441         586.0  74.8        57.7

What this tells us. The rush hours of working days (7:00 to 9:00 and 16:00 to 19:00) have the largest errors: 74.8 bikes per hour on average. There the model is too low by 57.7 bikes. These are only 441 of the 2,208 hours, but they are exactly the hours when stations run empty. As a share of the demand, the error is about 13% in rush hours and 17% at night. So the model is not worse at peaks in relative terms, but its misses there are large in bikes, and they point the wrong way.

In the notebook, you make the same table by weather code. Predict first: how much larger is the error in rain?

Show the answer

weather_slices = slice_table(valid, valid_forecast, valid["weathersit"])
print(weather_slices.to_string())

The output:

            hours  mean_rentals   MAE  mean_error
weathersit
1            1594         295.5  34.3        20.9
2             484         279.6  41.3        16.9
3             130         234.0  82.9        48.4

What this tells us. In light rain or snow (code 3), the MAE is 82.9: more than twice the 34.3 of clear weather, with only 130 hours. Rain is hard to place within an hour, and a shower at 17:00 changes the evening commute. In the real system, the forecast weather will add its own error on top. No hour in valid has code 4.

Last, the worst days of valid, by mean error per hour:

days = valid.assign(abs_error=np.abs(valid[TARGET] - valid_forecast), forecast=valid_forecast)
worst = days.groupby("dteday").agg(rentals=(TARGET, "sum"), forecast=("forecast", "sum"), MAE=("abs_error", "mean"))
print(worst.sort_values("MAE", ascending=False).head(4).round(0).to_string())

The output:

            rentals  forecast   MAE
dteday
2012-07-04     7403    5948.0  74.0
2012-09-08     5976    5370.0  68.0
2012-08-26     5255    4476.0  60.0
2012-07-21     4459    4304.0  58.0

What this tells us. The worst day is 4 July 2012, Independence Day: 7,403 rentals against a forecast of 5,948. The model only knows “holiday: yes”. It does not know that this holiday has fireworks in the city centre, and thousands of people ride there. The other bad days are weekends with fine weather. A list of special events would help more than a new algorithm.

The final test, once

The test quarter has been sealed since the split. Daniel writes down the plan before he runs anything. He refits E09 on train and valid together, so that the model knows the most recent level. Then he forecasts October to December 2012 once, and reports the result whatever it is. No setting changes after this point. You saw why in Tune hyperparameters.

train_and_valid = pd.concat([train, valid])
final = boosting_pipeline(chosen).fit(train_and_valid, train_and_valid[TARGET])
test_forecast = final.predict(test)
final_result = report(test[TARGET], test_forecast)
print("final model on test:", final_result)

The output:

final model on test: {'mean error': 15.3, 'MAE': 44.5, 'RMSE': 69.6}

For comparison, the two baselines, made in the same way from the data before October:

recent = train_and_valid[train_and_valid["dteday"] >= "2012-08-06"]  # the last 8 weeks before test
profile = recent.groupby(["workingday", "hr"])[TARGET].mean()
test_profile = profile.reindex(pd.MultiIndex.from_frame(test[["workingday", "hr"]])).to_numpy()
print("mean of train and valid:", report(test[TARGET], np.full(len(test), train_and_valid[TARGET].mean())))
print("hour profile:           ", report(test[TARGET], test_profile))

The output:

mean of train and valid: {'mean error': 34.0, 'MAE': 156.9, 'RMSE': 204.4}
hour profile:            {'mean error': -78.1, 'MAE': 94.1, 'RMSE': 150.0}

What this tells us. On the test quarter, the model is off by 44.5 bikes per hour, and too low by 15.3 on average. That is worse than on valid (38.7), which is normal. Autumn is a new season for the model, and nobody used the test set to choose anything. The hour profile does much worse in the autumn (94.1), because its summer shape is wrong when the days get colder and shorter. The model halves the error of the best simple plan.

The chart shows the week of Hurricane Sandy, hour by hour.

Real rentals
Forecast (dashed)

Real and forecast rentals for each hour from Friday 26 October to Thursday 1 November 2012, the week of Hurricane Sandy. The system had no rentals for most of 29 and 30 October.

Text description of the chart

A line chart of one week of hours, from Friday 26 October to Thursday 1 November 2012, with rentals per hour from 0 to 900. On Friday, Wednesday and Thursday, both lines show the two commuter peaks of a working day. The real peaks are higher than the forecast: at 17:00 on Friday, 817 real against 607 forecast. Saturday has one wide hill, with the forecast close to the real values. On Sunday 28 October, the storm arrives and both lines are lower. Monday 29 October has only one hour with data, 22 rentals. On Tuesday 30 October, the system reopens for 11 hours with 1,096 rentals, while the forecast for those hours adds up to 1,725. The lines jump over the hours that have no data.

What this tells us. On normal days, the forecast follows the shape of the day, but it is too low at the commuter peaks. After the storm, the system reopened on Tuesday for 11 hours. The forecast expected 1,725 rentals in those hours, and only 1,096 came. The model has never seen a city that is just reopening.

The slices of the test quarter are a report, not a reason to change the model:

print(slice_table(test, test_forecast, time_of_day(test)).to_string())
sandy = test["dteday"].isin(pd.to_datetime(["2012-10-29", "2012-10-30"]))
print("Hurricane Sandy, 29-30 Oct:", int(sandy.sum()), "hours, MAE",
      round(float(np.abs(test[TARGET] - test_forecast)[sandy].mean()), 1))
by_month = slice_table(test, test_forecast, test["dteday"].dt.month)
print(by_month[["hours", "mean_rentals", "MAE", "mean_error"]].to_string())

The output:

                hours  mean_rentals   MAE  mean_error
night (20-6)      991          80.8  21.5         9.1
other day hour    753         267.7  49.9         9.6
rush hour         424         456.7  88.7        39.9
Hurricane Sandy, 29-30 Oct: 12 hours, MAE 79.8
        hours  mean_rentals   MAE  mean_error
dteday
10        708         280.8  40.1        15.9
11        718         212.6  44.2        21.3
12        742         166.7  48.9         8.8

What this tells us. Rush hours are again the weakest slice (MAE 88.7, too low by 39.9). The 12 hours of Sandy have an MAE of 79.8: the model knew the rain but not that the city was shut. Each month of the quarter is too low on average, which is the growth of the system again.

Does the plan factor of 1.26 still cover 9 of 10 busy hours?

busy = test_forecast >= 100
covered = (test[TARGET].to_numpy()[busy] <= plan_factor * test_forecast[busy]).mean()
print("busy test hours:", int(busy.sum()), "| covered by forecast x", plan_factor, ":", round(float(covered), 3))

The output:

busy test hours: 1365 | covered by forecast x 1.26 : 0.788

What this tells us. It covers only 78.8% of the busy test hours, not 89.1%. The factor was chosen in the summer, and the autumn errors are larger. This is the most useful finding of the test: a safety margin is not a constant. It must be checked every month on the latest weeks, and changed when the coverage falls.

Last, save the final pipeline and look at the whole log, as in Package artifacts:

import joblib

joblib.dump(final, "bike_demand_model.joblib")
show_log()

The output:

 id                                                      change  valid_bias  valid_mae  valid_rmse decision
E01                  Baseline: the mean of train for every hour       120.6      197.8       260.7 baseline
E02 Baseline: mean per hour and day type, last 8 weeks of train        14.8       52.2        79.8 baseline
E03                Ridge, one-hot hour, weekday, month, weather        17.5      111.4       146.1   reject
E04                         Poisson regression, the same inputs       -28.0       93.0       131.9   reject
E05                        E04 + hour_type (hour x working day)       -27.9       54.8        87.3   reject
E06                  Gradient boosting, squared error, defaults        20.1       40.9        60.8   reject
E07                   Gradient boosting, Poisson loss, defaults        21.4       40.3        62.4     keep
E08          E07 + rentals in the same hour 2 and 7 days before        19.5       45.3        72.7   reject
E09                   E07 with learning_rate 0.05 and 400 trees        21.6       38.7        61.0     keep

What this tells us. The log has nine rows and the reason for every decision. Together with the saved pipeline, it lets another analyst repeat the work, or retrain next month in the same way.

Recommendations

Daniel writes a one-page note for Amara, in her words, not in model words:

Use the forecast to plan the busy hours of working days. It is off by about 45 bikes per hour in autumn, half the error of the 8-week profile you use today.
Plan with forecast × 1.26 for any hour with a forecast of 100 or more. In the summer, that covered 9 of 10 busy hours. In the autumn it covered 8 of 10. So at the start of every month, we will check the factor again on the last 4 weeks.
Expect the forecast to be too low, by about 15 bikes per hour, because the system keeps growing. We retrain it every month with the latest data.
Add special events by hand. The model does not know about fireworks, races or marches. On those days, your experience beats the model.
On days with severe weather warnings, decide by hand. The model knows rain, not a hurricane or a closed city.

And what not to use the model for. It does not tell you which station will run empty, because the data has only the whole system. It is not a tool to plan how many bikes to buy in the next years: it has seen only two years of fast growth. It does not forecast more than one day ahead, and it must not be used for any decision about individual staff.

Limits and next steps

Every result above has limits. Write them down before anyone else finds them:

Weather forecast, not observed weather. The data has the weather that happened. With a real forecast made the evening before, the errors will be larger, mostly on rainy days.
The whole system, not stations. Rebalancing needs a forecast per station, or per area. That needs the trip records of each station, which this dataset does not have.
Two years of fast growth. The model learns a level that keeps changing. Monthly retraining and a monthly check of the bias and of the plan factor are part of the model, not extras.
Missing hours. 165 hours have no row, most of them closures. The model never learned what a closure looks like, so it cannot forecast one.
Odd sensor values. Wind speed is 0 in 2,180 hours and humidity in 22. A next step is to replace them with values from another weather source.

The model card, as in Document use boundaries, in short:

Section	Content
Intended use	Hourly rentals of the whole system for the next day, made at 20:00, for staff and van planning.
Users	The operations team.
Not for	Station-level rebalancing, long-term planning, decisions about individual staff, severe weather or city closures.
Data	UCI Bike Sharing, Washington, D.C., 2011 and 2012, CC BY 4.0. Train to June 2012, valid July to September 2012, test October to December 2012.
Model	Gradient boosting with the Poisson loss, learning rate 0.05, 400 trees, 11 calendar and weather inputs.
Evaluation	Valid MAE 38.7. Test (once): MAE 44.5, RMSE 69.6, too low by 15.3. Hour-profile baseline on test: MAE 94.1.
Slices	Rush hours: test MAE 88.7, too low by 39.9. Light rain: about twice the clear-weather error on valid. Sandy: MAE 79.8.
Operating point	Plan factor 1.26 for forecasts of 100 or more; coverage 0.891 on valid, 0.788 on test.
Limits	Observed weather; whole system only; growth; holidays with events.
Owner and review	Daniel Kim; monthly review of the bias, the MAE of the last 4 weeks and the coverage of the plan factor. Retrain when the coverage of the last 4 weeks falls below 0.85.

Make it a portfolio piece

This case study is a good first project for a public portfolio: real data with an open licence, a real decision and an honest result. To make it yours, change at least one thing, for example a different split date, a station-level question, or another city’s open data. Then publish a small repository.

What to put in the repository:

bike-demand-forecast/
  README.md                problem, data and licence, approach, results, limits
  MODEL_CARD.md            the table above, in your words
  notebook.ipynb           the notebook, run from top to bottom, with outputs
  experiment_log.csv       every attempt, kept and rejected
  requirements.txt         pandas==3.0.6, scikit-learn==1.9.1, matplotlib

The README is the page that most people read. A README template:

# Hourly demand forecast for a bike-share system

## Problem
An operations team plans staff and vans each evening for the next day. Goal: forecast rentals per hour.

## Data and licence
UCI Bike Sharing (Fanaee-T, 2013), doi:10.24432/C5W894, CC BY 4.0. 17,379 hours, 2011-2012.

## Approach
Split by time (train to June 2012, valid Q3, test Q4). Removed the leak (casual + registered = target).
Baselines, then ridge, Poisson regression and gradient boosting with a Poisson loss. 9 experiments, logged.

## Results (test set, used once)
Gradient boosting: MAE 44.5 bikes per hour. Best simple plan (8-week hour profile): MAE 94.1.
A random split would have reported MAE 27: it is not a valid test for a forecast.

## Limits
Observed weather instead of forecast weather; whole system only; strong growth needs monthly retraining.

How to describe it in an interview, in two sentences: “I forecast hourly bike-share demand for next-day staff planning. A random split made the error look 30% smaller than it was, so I tested by time; the model halved the planner’s error, and the test showed that a summer safety margin fails in autumn.”

Go further:

Use a real weather forecast instead of the observed weather, for example from a public weather archive, and measure how much the error grows.
Forecast each station with the public trip data of a bike-share system, and plan the vans per area.
Replace the plan factor with a model that forecasts the 90th percentile directly: HistGradientBoostingRegressor(loss="quantile", quantile=0.9).
Add a calendar of public events and check the error on those days.
Run a monthly backtest over 2012 and plot the MAE and the bias by month.

Knowledge check

Check your understanding

Lena's random split gives an MAE of 27 bikes per hour, and Daniel's time split gives 38.7. Both use the same model. Explain the difference, and say which number Amara should expect.

Show the answer

With a random split, every test hour has its neighbours in train: the hour before and after, and the same weeks of growth and weather. The model learns the level of the very days it is tested on. Amara forecasts at 20:00 for tomorrow, when none of tomorrow’s hours are known. The time split copies that situation, so 38.7 on valid, and 44.5 on the later test quarter, are the numbers to expect. Even those are a little optimistic, because the data has the observed weather, not the forecast.

Check your understanding

Amara asks for the error as a percentage, because 'off by 45 bikes' means nothing to the board. What goes wrong with MAPE on this data, and what do you report instead?

Show the answer

MAPE divides each miss by the real count. At night, an hour can have 1 or 2 rentals, so a miss of 3 bikes becomes 150% or 300%. These hours then dominate the mean. MAPE then says the forecast is terrible when it is fine where it matters. Instead, report the MAE as a share of the mean demand: 44.5 out of 219.2 rentals per hour, about 20% on test. Or report MAPE only for the busy hours, and say which one you used.

Check your understanding

On the test quarter, the plan factor of 1.26 covered only 79% of busy hours instead of 90%. Daniel wants to choose a new factor on the test quarter. What do you advise?

Show the answer

Do not choose it on the test quarter. That would turn the test into a second validation set, and the 79% would no longer be an honest result. The finding is that the factor changes with the season. Write a rule. At the start of each month, choose the factor that covers 9 of 10 busy hours in the last 4 weeks of real data. At the end of the month, check its coverage. In live use, the newest weeks play the role of the validation set.

Recap

The problem was framed from the decision: one hour of the whole system, forecast at 20:00 the evening before, so casual and registered are leaks.
The split followed time, because the model forecasts the future. A random split reported MAE 27, a number Amara would never see.
The hour profile of the last 8 weeks was a strong baseline (valid MAE 52.2). Linear models did not beat it until they learned that counts multiply and that each kind of day has its own shape.
Gradient boosting with a Poisson loss won (valid MAE 38.7, never below zero). Past rentals made it worse, and the log keeps that result.
A plan factor of 1.26 turned the forecast into a plan. On the test quarter, used once, the model was off by 44.5 bikes per hour against 94.1 for the profile. But the factor covered only 79% of busy hours, so check margins every month.

Next, a second regression case study: a lender’s first estimate of a house price. The prices are skewed, errors are percentages, and the market moves under the model.

Next: Estimate house prices