Prévoir la demande de vélopartage

Étude de cas 1
120 min
Facultative · 1 étude de cas requise pour le certificat

Étude de cas 1 sur 7 · Apprentissage automatique : du problème à un modèle fiablePrévoir la demande de vélopartage

Aperçu du cours

1. Comprendre l'apprentissage automatique

2. Cadrer le problème de prédiction

3. Évaluer les modèles supervisés

4. Créer les premiers modèles

5. Essayer d'autres familles de modèles

6. Régler et finaliser le modèle

7. Comprendre et choisir le modèle

8. Préparer la passation du modèle

9. Études de cas

Ce que vous saurez faire

Partez d’une vraie question d’affaires pour arriver à une prévision testée et à une recommandation. Vous prévoyez combien de vélos un réseau urbain de vélopartage loue chaque heure. Vous voyez aussi pourquoi la réponse honnête est moins belle que la première.

Étude de casPièce de portfolio

Jeu de données: Bike Sharing, Washington DC 2011–2012 (UCI)
Licence: CC BY 4.0
Durée: Environ 2 h
S’exécute dans: Colab, Kaggle ou VS Code
Compétences utilisées: Détection des fuites
Division par temps
Modèles de Poisson
Erreurs par sous-groupe

Ouvrir dans Colab Ouvrir dans Kaggle Voir sur GitHub

Les carnets sont en anglais. Colab et Kaggle ouvrent votre propre copie. Pour utiliser VS Code, téléchargez le carnet depuis GitHub. Vous êtes bloqué ? Consultez le carnet de solution.

Le carnet (notebook) Jupyter contient toutes les cellules de code de cette page, trois tâches avec vérifications et les consignes. Ses cellules de préparation installent les paquets et téléchargent les données. Aucune configuration locale n’est donc nécessaire. Son exécution prend environ une minute. Lisez le récit ici et exécutez le code dans le carnet.

Avant de commencerCe qu’il faut savoir et avoir sous la mainEnviron 120 minutes

Ce que vous devez savoir

Cette étude de cas utilise tout le cours. Voici les leçons sur lesquelles elle s’appuie le plus :

Pourquoi évaluer une prévision sur des données ultérieures et utiliser les données de test une seule fois, dans Données d’entraînement, de validation et de test : comment séparer : au hasard ou selon le temps.
L’erreur moyenne, la MAE et la RMSE, et pourquoi la MAPE ne fonctionne plus près de zéro, dans Les métriques de régression.
Une référence de base qu’un modèle doit dépasser, dans Définir une référence de base.
Les pipelines avec encodage un à un et mise à l’échelle, dans Prétraiter avec un pipeline.
Le boosting (amélioration successive) par gradient, dans Les méthodes d’ensemble : le boosting : apprendre des erreurs, et le journal d’expériences, dans L’apprentissage automatique est itératif.
Les fuites, dans Sélectionner les caractéristiques utiles, et les sous-groupes avec leurs effectifs, dans L’analyse par sous-groupes.

Ce qu’il vous faut sur votre ordinateur (seulement pour l’exécution locale; le carnet n’en a pas besoin)

Un environnement Python avec pandas 3.0.6, scikit-learn 1.9.1 et matplotlib, comme le .venv de votre projet dans Définir une référence de base : préparer votre projet. Téléchargez le carnet depuis l’encadré ci-dessus dans un dossier de votre choix. Ouvrez-le comme dans Définir une référence de base : ouvrir un carnet dans VS Code. Le carnet télécharge lui-même les données dans un nouveau dossier bike-demand/. Pour recommencer, supprimez ce dossier.

Testé sur macOS avec Python 3.14.6, pandas 3.0.6, scikit-learn 1.9.1 et matplotlib 3.11.2. Le carnet a été testé dans Jupyter. Windows n’a pas été testé.

La situation

Amara Okafor dirige les opérations d’un réseau urbain de vélopartage. Il compte quelques centaines de stations où les gens prennent un vélo et le laissent à une autre station. Chaque soir, elle planifie le lendemain. Elle décide combien d’employés travaillent aux stations les plus fréquentées. Elle décide aussi combien de fourgonnettes déplacent les vélos des stations pleines vers les stations vides.

Les deux types d’erreurs lui coûtent cher. Si elle prévoit trop peu de ressources pour une heure achalandée, les gens trouvent des stations vides à 8 h. Certains cessent d’utiliser le réseau. Si elle prévoit trop de ressources, elle paie des conducteurs et des fourgonnettes qui attendent. La première erreur est plus grave : perdre un client coûte plus cher qu’une heure de fourgonnette. Aujourd’hui, elle planifie selon son expérience et les chiffres de la semaine précédente.

Elle pose une question précise à Daniel Kim, l’analyste de données : « Combien de vélos les gens loueront-ils chaque heure demain? » Elle a besoin de la réponse à 20 h, avec les prévisions météo du lendemain. Elle peut payer environ deux fourgonnettes supplémentaires aux heures de pointe. Elle veut donc savoir quand les pointes surviendront et quelle sera leur ampleur.

Ce que vous allez créer

Vous allez créer une prévision de la demande horaire du lendemain et une règle qui la transforme en plan. Vous la comparerez honnêtement à ce qu’une personne peut planifier sans apprentissage automatique, et vous noterez ses échecs. À la fin, vous aurez les éléments d’un projet pour votre portfolio : un carnet, un journal de neuf expériences, un modèle sauvegardé, une courte fiche de modèle et une recommandation.

L’étude de cas utilise les compétences de tout le cours. Vous cadrez le problème et détectez une fuite (modules 2 et 6). Vous divisez les données par temps et choisissez des métriques de régression (module 3). Vous créez des pipelines (module 4), utilisez le boosting et tenez un journal d’expériences (module 5). Vous réglez un paramètre (module 6). Ensuite, vous vérifiez les sous-groupes (module 7) et rédigez une fiche de modèle (module 8). Deux idées sont nouvelles, et la page les explique au moment voulu : la perte de Poisson pour les dénombrements et un facteur de planification qui transforme une prévision en décision.

Les données

Les données proviennent du jeu de données Bike Sharing, dans le UCI Machine Learning Repository. Il contient chaque heure du réseau Capital Bikeshare à Washington, D.C., en 2011 et 2012, avec la météo correspondante. Hadi Fanaee-T et João Gama l’ont créé à partir des registres publics de trajets du réseau. Ils ont ajouté la météo de freemeteo.com et la liste des jours fériés de la ville.

Source : Fanaee-T, H. (2013). Bike Sharing [jeu de données]. UCI Machine Learning Repository. doi:10.24432/C5W894. Sous licence CC BY 4.0, qui permet la réutilisation avec attribution. Les auteurs vous demandent de citer leur article : Fanaee-T, H. et Gama, J. (2013), Event labeling combining ensemble detectors and background knowledge, Progress in Artificial Intelligence, doi:10.1007/s13748-013-0040-3. Nextia héberge le fichier hour.csv inchangé et comprimé, pour que le carnet reçoive les mêmes octets à chaque téléchargement.

Le réseau d’Amara n’est pas Capital Bikeshare. Amara et Daniel sont des personnages inventés pour relier les données à une décision. Les chiffres sont réels.

Charger les données et regarder trois lignes

Chaque bloc de code de cette page est une cellule du carnet. Les cellules de préparation ont déjà chargé pandas sous pd et numpy sous np, et téléchargé le fichier dans le dossier de travail.

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())

Le résultat :

(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

Ce que cela nous apprend. Le fichier contient 17 379 lignes et 17 colonnes. Une ligne représente une heure pour tout le réseau : à minuit le 1er janvier 2011, il y avait 16 locations. La cible est cnt, le nombre de locations dans l’heure. casual (occasionnels) et registered (abonnés) répartissent ce nombre entre utilisateurs sans abonnement et avec abonnement. Regardez temp (température) : 0,24 n’est pas une température en degrés. Les auteurs ont divisé chaque température par 41, donc 0,24 signifie environ 10 °C.

Les colonnes utilisées

Avant de créer un modèle, notez ce que signifie chaque colonne. Voici les colonnes utilisées dans cette étude de cas, en mots simples :

Colonne	Signification
`dteday`, `hr`	la date et l’heure (0 à 23)
`weekday`, `workingday`, `holiday`	jour de la semaine (0 = dimanche); 1 si ce n’est ni une fin de semaine ni un jour férié; 1 un jour férié
`mnth`, `yr`	mois (1 à 12); année (0 = 2011, 1 = 2012)
`weathersit`	1 ciel dégagé ou partiellement nuageux; 2 brume ou ciel nuageux; 3 pluie ou neige légère; 4 forte pluie, glace ou neige
`temp`, `atemp`	température et température « ressentie », mises à l’échelle (× 41 et × 50 donnent des °C)
`hum`, `windspeed`	humidité (× 100 donne des %) et vitesse du vent, mises à l’échelle
`cnt`	la cible : locations dans l’heure
`casual`, `registered`	locations par les utilisateurs occasionnels et par les abonnés

Nous écartons deux colonnes. instant est seulement le numéro de ligne, et season répète ce que mnth indique déjà.

Premières vérifications

Ensuite, faites les vérifications nécessaires pour tout nouveau jeu de données : cellules vides, lignes en double, périodes manquantes et codes suspects.

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()))

Le résultat :

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

Ce que cela nous apprend. Aucune cellule n’est vide et aucune heure n’apparaît deux fois. Mais deux années comptent 17 544 heures, et le fichier contient 17 379 lignes. Donc 165 heures n’ont aucune ligne. Le code météo 4 apparaît seulement à 3 heures, trop peu pour qu’un modèle puisse en apprendre quelque chose. Une humidité de 0 est impossible à Washington : ces 22 valeurs correspondent donc à des mesures manquantes. Une vitesse du vent exactement égale à 0 pendant 2 180 heures est aussi suspecte. Elle signifie probablement « non mesurée ». Nous gardons ces lignes et retenons que la vitesse du vent n’est pas fiable.

Où sont les heures manquantes? Comptez-les par jour :

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())

Le résultat :

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

Ce que cela nous apprend. Les plus grandes périodes manquantes sont les 29 et 30 octobre 2012. C’était l’ouragan Sandy : le réseau avait fermé. Les autres périodes manquantes, en janvier 2011, correspondent à des journées de neige et à des problèmes de données. Les données n’ont donc aucune ligne pour une heure où personne ne pouvait louer de vélo. Ces deux journées de Sandy se trouvent au dernier trimestre de 2012, qui sera notre période de test.

La cible et une fuite

Regardez maintenant la cible et vérifiez un soupçon : est-ce que casual et registered donnent cnt quand on les additionne?

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

Le résultat :

{'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

Ce que cela nous apprend. Une heure typique compte environ 142 locations, et l’heure la plus achalandée en compte 977. Le plus petit nombre est 1, jamais 0 : une heure sans location n’a probablement aucune ligne. Et casual + registered est égal à la cible dans chaque ligne. Ces deux colonnes sont comptées dans la même heure que la cible. À 20 h, leurs valeurs du lendemain sont donc inconnues. Elles constituent une fuite : une information indisponible au moment de la prédiction, qui fait paraître un modèle beaucoup meilleur qu’il ne peut l’être. Vous avez vu les fuites dans Sélectionner les caractéristiques utiles. Ces colonnes ne font jamais partie des entrées.

Cadrer le problème

Décrivez maintenant le problème selon les termes de Traduire un besoin, avant de créer un modèle :

Décision : combien d’employés et de fourgonnettes Amara prévoit pour chaque heure du lendemain.
Unité de prédiction : une heure pour tout le réseau.
Cible : cnt, les locations dans cette heure. C’est un nombre, donc il s’agit de régression.
Moment de la prédiction : 20 h la veille. Tout ce qu’une entrée utilise doit être connu à ce moment.
Entrées disponibles à ce moment : le calendrier (heure, jour de la semaine, mois, année, jour férié, jour ouvrable) et les prévisions météo du lendemain.
Informations indisponibles : casual et registered (la fuite), et les locations des heures précédentes du lendemain.
Personnes : aucune colonne ne décrit une personne, donc aucun sous-groupe ne permet d’analyser l’équité entre groupes. Les sous-groupes importants sont les heures, la météo et les journées particulières.

Une limite apparaît déjà. Les données contiennent la météo réelle, pas les prévisions qu’Amara aura à 20 h. Les prévisions météo du lendemain sont généralement proches de la réalité, mais elles peuvent manquer une tempête. Chaque erreur sur cette page est donc un peu plus faible que celle qu’Amara observera. Nous le notons maintenant, puis de nouveau dans les limites.

Explorer ce dont la décision a besoin

Nous n’explorons pas toutes les colonnes. Nous posons trois questions dont dépend le plan d’Amara. À quoi ressemble une journée? La demande change-t-elle? Quelle importance a la météo?

Le profil d’une journée

Amara planifie heure par heure, donc la première question porte sur le profil d’une journée. Calculez la moyenne des locations pour chaque heure, séparément pour les jours ouvrables et les autres jours :

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())

Le résultat (cinq des 24 heures) :

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

Le graphique montre les 24 heures.

Jour ouvrable
Fin de semaine ou jour férié

Locations moyennes pour chaque heure de la journée, en 2011 et 2012. Les jours ouvrables ont deux pointes liées aux déplacements domicile-travail; les autres jours ont une large pointe vers midi.

Description textuelle du graphique

Un graphique en courbes. L’axe horizontal représente l’heure de la journée, de 0 à 23. L’axe vertical représente les locations moyennes par heure, de 0 à 600. La courbe des jours ouvrables est basse la nuit, avec environ 5 locations à 3 h et à 4 h. Elle monte rapidement jusqu’à une pointe de 477 à 8 h, puis redescend entre environ 135 et 200 au milieu de la journée. Elle remonte ensuite à son maximum, 525 à 17 h et 492 à 18 h, puis redescend à 89 à 23 h. La courbe des fins de semaine et des jours fériés reste basse jusqu’à 7 h. Elle monte jusqu’à une large pointe d’environ 355 à 373 entre 12 h et 16 h, puis redescend à 86 à 23 h. À 8 h, un jour ouvrable compte environ quatre fois et demie les locations d’un autre jour. À 13 h, les autres jours en comptent presque deux fois plus.

Ce que cela nous apprend. Un jour ouvrable a deux pointes, à 8 h et entre 17 h et 18 h : les gens vont travailler et reviennent. Les fins de semaine et les jours fériés ont une large pointe vers midi : les gens roulent pour le plaisir. L’heure seule ne suffit donc pas. Une prévision doit connaître le type de journée. Elle doit aussi bien prévoir les deux pointes des déplacements domicile-travail, car c’est à ces heures que les stations se vident.

La tendance

Ensuite, la demande change-t-elle avec le temps? Comparez les locations moyennes par heure des deux années :

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()))

Le résultat :

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

La moyenne est passée de 143,8 à 234,7 locations par heure, soit une hausse de 63 % en un an. Le graphique montre les locations moyennes par jour pour chaque semaine. Les couleurs correspondent aux parties de la division présentée dans la prochaine section.

Entraînement
Validation
Test

Locations moyennes par jour, semaine par semaine, de janvier 2011 à décembre 2012. Chaque été est plus achalandé que chaque hiver, et 2012 est plus achalandée que 2011 à toutes les saisons.

Description textuelle du graphique

Un graphique en courbes de 106 semaines, de la fin décembre 2010 à la fin décembre 2012. L’axe vertical représente les locations moyennes par jour, de 0 à 8 000. La partie d’entraînement, jusqu’en juin 2012, commence sous 1 000 locations par jour en janvier 2011. Elle monte entre environ 4 000 et 5 000 à l’été 2011, puis redescend à environ 3 000 en hiver. Elle remonte ensuite à environ 7 000 en juin 2012. La partie de validation, de juillet à septembre 2012, est la plus élevée : entre 6 036 et 7 860 locations par jour. La partie de test, d’octobre à décembre 2012, commence entre environ 6 500 et 7 000. Elle chute à 4 109 pendant la semaine de l’ouragan Sandy. Elle reste entre environ 3 900 et 5 800 en novembre et au début décembre. Elle redescend à 1 531 pendant la semaine de Noël.

Ce que cela nous apprend. Deux tendances se superposent : les saisons et une forte croissance. L’été 2012 est beaucoup plus achalandé que l’été 2011, car le réseau a ajouté des stations et des abonnés. C’est important pour tous les modèles. Un modèle apprend le niveau du passé, donc ses prévisions auront tendance à être trop basses pour l’avenir. Et le dernier trimestre de 2012 combine une année achalandée et une saison où la demande baisse.

La météo

Enfin, à quel point la météo change-t-elle la demande?

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

Le résultat :

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

Ce que cela nous apprend. Sous une pluie ou une neige légère (code 3), le réseau loue environ 112 vélos par heure. C’est presque la moitié des 205 locations par temps dégagé. Les prévisions météo du lendemain sont donc une entrée importante. Les heures pluvieuses peuvent être difficiles à prévoir, car la pluie commence et s’arrête souvent au cours d’une heure.

Division et référence de base

Une division chronologique

Dans Données d’entraînement, de validation et de test, vous divisez les données chronologiquement quand le modèle doit prédire l’avenir. C’est exactement le cas d’Amara : chaque soir, le modèle connaît le passé et prévoit demain. La division suit donc le calendrier :

Entraînement : de janvier 2011 à la fin de juin 2012. Les modèles apprennent à partir de ces données.
Validation : de juillet à septembre 2012. Nous les utilisons pour comparer les modèles et prendre toutes les décisions.
Test : d’octobre à décembre 2012. Ces données restent fermées jusqu’au test final.

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))

La sortie :

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

Ce que cela nous dit. Les données d’entraînement comptent 13 003 heures (18 mois). Les données de validation et de test en comptent environ 2 200 chacune (3 mois). La validation est beaucoup plus achalandée que l’entraînement : 288,4 locations par heure contre 167,7. Cela combine l’été et la croissance. Un modèle entraîné surtout sur des mois plus calmes sera testé sur le trimestre le plus achalandé.

Comment nous mesurons les erreurs

Quelle métrique convient à la décision d’Amara? Elle interprète les erreurs en nombre de vélos. Le chiffre principal est donc la MAE, l’écart typique en vélos par heure. Une prévision trop basse laisse les stations vides, donc le sens de l’erreur compte. Nous indiquons toujours aussi l’erreur moyenne : une valeur positive signifie que la prévision est trop basse. Les grands écarts aux pointes sont les pires. Nous indiquons donc aussi la RMSE, qui les pénalise. Nous n’utilisons pas la MAPE : la nuit, une heure peut compter 2 locations, et la MAPE fonctionne mal près de zéro. Vous l’avez vu dans Les métriques de régression.

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)}

Deux références de base

Créer un modèle vaut la peine seulement s’il fait mieux qu’Amara sans modèle. Vous l’avez appris dans Définir une référence de base. Daniel essaie deux références de base. La première utilise la moyenne des données d’entraînement pour chaque heure : c’est le minimum à dépasser. La seconde reproduit ce qu’une personne attentive ferait dans un tableur : le profil horaire. C’est la moyenne des locations pour chaque heure des jours de travail et des autres jours. Elle couvre les 8 dernières semaines d’entraînement. Les semaines récentes suivent mieux la croissance que l’historique complet.

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)

La sortie :

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}

Ce que cela nous dit. La moyenne des données d’entraînement se trompe d’environ 198 vélos pour une heure typique. Elle ignore le profil de la journée. Le profil horaire se trompe d’environ 52 vélos par heure. Il est trop bas d’environ 15 vélos en moyenne, car la croissance estivale s’est poursuivie après les 8 semaines. Le profil est une référence de base solide, et il ne coûte rien. Un modèle doit faire mieux que 52,2, pas que 197,8.

Le journal des expériences

À partir d’ici, chaque essai entre dans un journal des expériences, comme dans L’apprentissage automatique est itératif. Les colonnes sont les mêmes, avec les métriques de ce problème : id (identifiant), date, change (changement), valid_bias (erreur moyenne de validation), valid_mae (MAE de validation), valid_rmse (RMSE de validation), decision (décision) et reason (raison). Avec des données chronologiques, le trimestre de validation joue le rôle de la validation croisée utilisée pour les billets de soutien.

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()

La sortie :

 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

Ce que cela nous dit. Le journal contient maintenant les deux références de base, avec la raison de chacune. Chaque modèle doit désormais battre E02, avec une MAE de validation inférieure à 52,2.

Modèles

Daniel essaie trois types de modèles, du plus simple au plus complexe, et conserve chaque résultat. Chaque modèle est un pipeline scikit-learn : le prétraitement et le modèle dans un seul objet, ajusté uniquement sur les données d’entraînement.

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)])

Les deux fonctions choisissent leurs propres colonnes. linear_pipeline() encode l’heure, le jour de la semaine, le mois et le code météo par encodage disjonctif. Elle met aussi les nombres à l’échelle, comme dans Prétraiter avec un pipeline. boosting_pipeline() conserve seulement les colonnes indiquées, car les arbres n’ont pas besoin de mise à l’échelle. Chaque pipeline sélectionne ses propres entrées. Vous pouvez donc lui passer le DataFrame complet, sans que casual ou registered s’y glisse par erreur.

Un modèle linéaire, et pourquoi il échoue

Le premier modèle est la régression ridge, une régression linéaire qui pénalise légèrement les grands coefficients. Vous avez découvert cette pénalité comme régularisation dans Contrôler la complexité.

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)

La sortie :

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

Ce que cela nous dit. La MAE est de 111,4 : deux fois pire que le profil horaire. Pourquoi? Comparez la prévision à la moyenne réelle pour quatre heures d’un jour de travail :

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())

La sortie :

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

Ce que cela nous dit. À 3 h, environ 7 personnes louent un vélo, et le modèle en prévoit 109. À 17 h, 813 personnes louent un vélo, et il en prévoit 518. Le modèle est trop uniforme : trop haut la nuit, trop bas aux pointes. La raison tient au fonctionnement d’un modèle linéaire. Il additionne les effets : « un jour de travail ajoute tant de vélos, la pluie en enlève tant, 8 h en ajoute tant ». Mais la demande ne fonctionne pas ainsi.

Les nombres se multiplient : le modèle de Poisson

Pensez à un comptoir de crème glacée. Par temps chaud, il vend deux fois plus que d’habitude, le matin comme l’après-midi. Il ne vend pas « 200 de plus » chaque heure. À 7 h, quand il en vend habituellement 10, il en vend 20. À 15 h, quand il en vend habituellement 300, il en vend 600. La météo et la saison multiplient la demande de chaque heure. C’est pareil pour les vélos : la pluie réduit les locations d’environ la moitié, aux heures calmes comme aux heures achalandées.

Un modèle où les effets se multiplient est un modèle linéaire appliqué au logarithme du nombre attendu. Sur l’échelle logarithmique, « × 2 » devient « + log 2 » : l’addition redevient donc une multiplication. Le modèle de comptage qui fait cela est la régression de Poisson. Il porte le nom de Siméon Poisson, un mathématicien français qui étudiait le nombre d’événements. Il étudiait, par exemple, le nombre de lettres perdues par la poste chaque jour. Ses prévisions ne sont jamais négatives, ce qui convient à un comptage.

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)

La sortie :

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

Ce que cela nous dit. La MAE passe de 111,4 à 93,0. C’est mieux, mais encore loin des 52,2 du profil. La prévision est maintenant trop haute en moyenne (erreur moyenne de −28,0). Le modèle utilise encore un seul profil quotidien pour tous les jours. Il ne peut pas savoir que 8 h est une pointe le lundi et une heure calme le dimanche. Daniel lui donne directement cette information : une catégorie pour chaque heure de chaque type de jour, soit 48 valeurs au total :

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)

La sortie :

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

Ce que cela nous dit. La MAE descend à 54,8. Le profil de la journée était l’élément manquant. Une combinaison de deux colonnes dont l’effet change lorsqu’elles agissent ensemble s’appelle une interaction. Les modèles linéaires exigent que vous la construisiez; nous verrons que les arbres la trouvent seuls. Malgré cela, le meilleur modèle linéaire reste un peu moins bon que le profil horaire. Ajoutons ces résultats au journal :

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()

La sortie :

 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

Ce que cela nous dit. Trois modèles linéaires, et aucun ne bat le profil horaire. Le journal conserve les trois, avec leur raison, pour éviter de les réessayer sans nouvelle idée.

Le boosting par gradient, avec la bonne perte

Ensuite vient le boosting par gradient, une méthode de renforcement présentée dans Les méthodes d’ensemble. Elle combine de nombreux petits arbres : chaque nouvel arbre corrige les erreurs des précédents. Un arbre peut d’abord diviser selon le type de jour, puis selon l’heure. Il trouve donc l’interaction sans aide. Le HistGradientBoostingRegressor de scikit-learn est la version pour la régression du classificateur que vous avez utilisé.

Le boosting mesure ses erreurs avec une perte, comme la régression linéaire dans Comprendre la régression. Par défaut, il utilise l’erreur quadratique. Daniel essaie aussi loss="poisson" : la même idée que la régression de Poisson. Cette perte évalue un écart par rapport au nombre prévu. Un écart de 20 vélos compte davantage dans une heure calme que dans une heure achalandée. Ses prévisions ne sont jamais négatives.

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()))

La sortie :

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

Ce que cela nous dit. Les deux modèles de boosting battent le profil horaire d’environ 12 vélos par heure : 40,9 et 40,3 contre 52,2. Le modèle à erreur quadratique prévoit moins de 0 vélo pendant 10 heures nocturnes. C’est impossible, et une personne qui planifie perd confiance dans l’outil en voyant cela. La perte de Poisson ne produit jamais ce résultat, et sa MAE est un peu plus faible. Daniel la conserve donc. Les deux modèles sont trop bas d’environ 20 vélos par heure en moyenne. La validation est plus achalandée que toutes leurs données d’entraînement.

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()

La sortie :

 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

Ce que cela nous dit. E07 est la première ligne marquée keep (conserver). C’est le premier modèle qui bat le profil horaire, de 12 vélos par heure.

À vous : les locations passées comme entrées

Daniel a une autre idée. Les locations récentes devraient nous renseigner sur demain. Quelles locations passées connaît-on à 20 h? Pas celles de la même heure la veille. À 20 h aujourd’hui, les heures de 20 h à minuit n’ont pas encore eu lieu. Pour prévoir demain à 21 h, il faudrait les locations d’aujourd’hui à 21 h. Les locations de la même heure 2 jours avant et 7 jours avant sont connues. Il ajoute donc les deux, avec une petite fonction auxiliaire :

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})

Dans le carnet (notebook) Jupyter, ajoutez les deux colonnes aux données d’entraînement et de validation. Ajustez le modèle de boosting de Poisson avec elles, puis consignez le résultat comme E08. Faites d’abord une prédiction : les locations passées seront-elles utiles?

Afficher la réponse

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.")

La sortie :

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

Ce que cela nous dit. Les locations passées rendent la prévision moins bonne : une MAE de 45,3 contre 40,3. Les locations d’il y a deux jours reflètent la météo de ce jour-là, un mardi pluvieux ou un samedi ensoleillé. Demain ne partage pas cette météo. Le modèle reçoit une entrée qui semble utile, mais qui est surtout du bruit pour cette décision. Une idée rejetée reste un résultat, donc elle reste dans le journal. Vous ne le sauriez pas sans essayer, et vous n’y feriez pas confiance sans le journal.

Régler un paramètre et le vérifier deux fois

La dernière expérience ralentit l’apprentissage : un taux d’apprentissage de 0,05 et 400 arbres, plutôt que 0,1 et au plus 100. Dans Régler les hyperparamètres, vous avez comparé les réglages avec cinq plis. Ici, les données sont chronologiques, et un trimestre de validation donne un seul résultat bruité. Daniel vérifie donc aussi le réglage sur un second trimestre. Il entraîne le modèle sur les données antérieures à avril 2012, puis prévoit d’avril à juin 2012. C’est un test rétrospectif : un test sur une période passée, comme si vous y étiez.

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)))

La sortie :

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}

Ce que cela nous dit. Le modèle plus lent gagne sur les deux trimestres. Il obtient 38,7 contre 40,3 en validation, et 37,6 contre 39,3 au printemps. Les gains sont petits, mais concordent : E09 devient donc le réglage final.

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()

La sortie :

 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

Le graphique résume le journal. E01, la moyenne des données d’entraînement, est omise : sa MAE de 197,8 écraserait toutes les autres barres.

MAE de validation de chaque expérience, de E02 à E09. Gris : le profil horaire de référence. Bleu : les modèles. Orange : le réglage final, E09.

Description textuelle du graphique

Un graphique à colonnes de huit expériences, avec la MAE de validation en vélos par heure, de 0 à 120. E02, le profil horaire, la référence de base : 52,2. E03, la régression ridge : 111,4, la barre la plus haute. E04, la régression de Poisson : 93,0. E05, la régression de Poisson avec le type d’heure : 54,8. E06, le boosting à erreur quadratique : 40,9. E07, le boosting avec la perte de Poisson : 40,3. E08, E07 avec les locations passées : 45,3. E09, le boosting plus lent : 38,7, la barre la plus basse. Seuls les quatre modèles de boosting sont sous la référence de base.

Ce que cela nous dit. Le journal raconte le travail en une seule image. Les modèles linéaires se sont approchés de la référence seulement après avoir appris le profil de chaque type de jour. Tous les modèles de boosting l’ont battue. Les deux réglages marqués keep sont les deux étapes qui ont aidé.

Erreur courante

La division aléatoire qui semble excellente

Symptôme : Lena, une nouvelle analyste de l’équipe, suit un tutoriel populaire. Elle mélange toutes les heures et en garde 20 % pour le test. Son modèle se trompe de 27 vélos par heure, bien mieux que les 38,7 de Daniel. Pour faire encore mieux, elle ajoute casual, et l’erreur descend à 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)))

La sortie :

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 : avec une division aléatoire, chaque heure de test a ses voisines dans les données d’entraînement. Il y a l’heure précédente, l’heure suivante, la même heure la veille. Le modèle apprend le niveau de ces mêmes journées, y compris la croissance de 2012 et chaque tempête. Amara n’a jamais cela : à 20 h, les heures voisines de demain n’existent pas encore. Avec casual, le modèle voit aussi une partie de la réponse. Les deux chiffres mesurent une situation qui ne se produira jamais.

Correction : divisez les données chronologiquement, selon l’utilisation prévue du modèle. Apprenez du passé et testez sur une période ultérieure. Retirez toute entrée inconnue au moment de la prédiction. Attendez-vous ensuite à une erreur plus grande, mais honnête : ici, 38,7 en validation, et davantage au test.

Conséquence dans le travail en IA : une prévision testée avec une division aléatoire semble excellente dans le rapport, puis déçoit dès le premier mois. L’équipe accuse alors le modèle, les données ou les utilisateurs, alors que la vraie cause était le test.

Métrique et point de fonctionnement

Le modèle choisi, E09, ajusté sur les données d’entraînement et vérifié sur les données de validation :

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))

La sortie :

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

Ce que cela nous dit. Ce sont les chiffres de E09 dans le journal. L’écart typique est de 38,7 vélos par heure, et la prévision est trop basse de 21,6 en moyenne.

Une prévision n’est pas encore un plan. Le modèle est trop bas en moyenne, et pour Amara, une pénurie est pire qu’un vélo de surplus. Elle ne planifie donc pas avec la prévision telle quelle. Pensez à une boulangère qui sait que son commerce vend environ 100 pains le samedi. Elle en prépare 120, car une rupture de stock fait perdre des clients, et quelques pains de surplus coûtent peu. La quantité supplémentaire est sa marge de sécurité.

Amara fait la même chose avec un facteur de planification. Aux heures achalandées, celles où la prévision atteint au moins 100 vélos, elle prévoit prévision × facteur. Daniel choisit le facteur sur les données de validation pour que le plan couvre 9 heures achalandées sur 10. C’est la version pour la régression du seuil dans Seuils et calibration. Ce nombre transforme la sortie du modèle en décision. On le choisit sur les données de validation selon le coût de chaque erreur. Un tel nombre est le point de fonctionnement.

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))

La sortie :

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

Ce que cela nous dit. Avec un facteur de 1,26, le plan couvre 89,1 % des 1 618 heures achalandées de validation. Pendant une heure achalandée de validation, cela représente environ 92 vélos de capacité supplémentaire en moyenne. C’est le prix de la marge de sécurité, et Amara peut le comparer aux fourgonnettes dont elle dispose.

Erreurs et sous-groupes

Où le modèle échoue-t-il? Comme dans L’analyse par sous-groupes, nous examinons des groupes d’heures de validation, avec leur nombre d’heures à côté de chaque erreur. Toutes les décisions utilisent la validation. Les sous-groupes de test viendront plus tard, uniquement comme rapport.

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())

La sortie :

                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

Ce que cela nous dit. Les heures de pointe des jours de travail (de 7 h à 9 h et de 16 h à 19 h) ont les plus grandes erreurs : 74,8 vélos par heure en moyenne. Le modèle y est trop bas de 57,7 vélos. Cela représente seulement 441 des 2 208 heures, mais ce sont exactement les heures où les stations se vident. Par rapport à la demande, l’erreur est d’environ 13 % aux heures de pointe et de 17 % la nuit. Le modèle n’est donc pas moins bon aux pointes en termes relatifs. Mais ses écarts y représentent beaucoup de vélos, et ils vont dans le mauvais sens.

Dans le carnet, faites le même tableau par code météo. Faites d’abord une prédiction : de combien l’erreur augmente-t-elle sous la pluie?

Afficher la réponse

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

La sortie :

            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

Ce que cela nous dit. Sous une pluie ou une neige légère (code 3), la MAE est de 82,9. C’est plus du double des 34,3 par temps clair, avec seulement 130 heures. Il est difficile de situer la pluie dans une heure, et une averse à 17 h change le trajet du soir. Dans le vrai système, les prévisions météo ajouteront leurs propres erreurs. Aucune heure de validation n’a le code 4.

Enfin, les pires journées de validation, selon l’erreur moyenne par heure :

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())

La sortie :

            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

Ce que cela nous dit. La pire journée est le 4 juillet 2012, jour de l’Indépendance : 7 403 locations contre une prévision de 5 948. Le modèle sait seulement « jour férié : oui ». Il ignore que cette fête comprend des feux d’artifice au centre-ville et que des milliers de personnes s’y rendent à vélo. Les autres mauvaises journées sont des fins de semaine de beau temps. Une liste d’événements spéciaux aiderait davantage qu’un nouvel algorithme.

Le test final, une seule fois

Le trimestre de test est resté fermé depuis la séparation des données. Daniel écrit le plan avant de lancer quoi que ce soit. Il réentraîne E09 sur les données d’entraînement et de validation réunies, pour que le modèle connaisse le niveau le plus récent. Ensuite, il prévoit les locations d’octobre à décembre 2012 une seule fois, et rapporte le résultat, quel qu’il soit. Aucun réglage ne change après ce point. Vous avez vu pourquoi dans Régler les hyperparamètres.

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)

Le résultat :

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

Pour comparer, voici les deux références de base, calculées de la même façon avec les données antérieures à octobre :

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))

Le résultat :

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}

Ce que cela nous apprend. Sur le trimestre de test, le modèle se trompe de 44,5 vélos par heure. Il sous-estime la demande de 15,3 vélos en moyenne. C’est pire que sur les données de validation (38,7), ce qui est normal. L’automne est une nouvelle saison pour le modèle, et personne n’a utilisé les données de test pour choisir quoi que ce soit. Le profil horaire fait beaucoup moins bien à l’automne (94,1). Sa forme estivale ne convient plus quand les journées deviennent plus froides et plus courtes. Le modèle réduit de moitié l’erreur du meilleur plan simple.

Le graphique montre la semaine de l’ouragan Sandy, heure par heure.

Locations réelles
Prévision (pointillés)

Locations réelles et prévues pour chaque heure du vendredi 26 octobre au jeudi 1 novembre 2012, la semaine de l'ouragan Sandy. Le réseau n'a enregistré aucune location pendant la majeure partie des 29 et 30 octobre.

Description textuelle du graphique

Un graphique en courbes couvre une semaine, heure par heure, du vendredi 26 octobre au jeudi 1 novembre 2012. L’axe des locations par heure va de 0 à 900. Le vendredi, le mercredi et le jeudi, les deux courbes montrent les deux pointes des déplacements domicile-travail d’une journée de travail. Les pointes réelles dépassent les prévisions : à 17 h le vendredi, 817 locations réelles contre 607 prévues. Le samedi présente une seule bosse large, et la prévision reste proche des valeurs réelles. Le dimanche 28 octobre, la tempête arrive et les deux courbes sont plus basses. Le lundi 29 octobre ne compte qu’une heure avec des données, soit 22 locations. Le mardi 30 octobre, le réseau rouvre pendant 11 heures et enregistre 1 096 locations. La prévision pour ces heures totalise 1 725 locations. Les courbes sautent les heures sans données.

Ce que cela nous apprend. Les jours normaux, la prévision suit la forme de la journée. Mais elle est trop basse aux pointes des déplacements domicile-travail. Après la tempête, le réseau a rouvert le mardi pendant 11 heures. La prévision annonçait 1 725 locations pendant ces heures, mais il n’y en a eu que 1 096. Le modèle n’a jamais vu une ville qui vient de reprendre ses activités.

Les sous-groupes du trimestre de test servent à produire un rapport, pas à modifier le modèle :

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())

Le résultat :

                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

Ce que cela nous apprend. Les heures de pointe restent le sous-groupe le plus faible (MAE de 88,7, sous-estimation de 39,9). Les 12 heures de Sandy ont une MAE de 79,8. Le modèle connaissait la pluie, mais pas l’arrêt des activités de la ville. Il sous-estime la demande moyenne de chaque mois du trimestre. C’est encore l’effet de la croissance du réseau.

Le facteur de planification de 1,26 couvre-t-il encore 9 heures achalandées sur 10?

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))

Le résultat :

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

Ce que cela nous apprend. Il ne couvre que 78,8 % des heures achalandées du test, et non 89,1 %. Le facteur a été choisi en été, et les erreurs sont plus grandes en automne. C’est la conclusion la plus utile du test : une marge de sécurité n’est pas une constante. Il faut la vérifier chaque mois sur les semaines les plus récentes, et la modifier quand la couverture baisse.

Enfin, enregistrez le pipeline final et examinez le journal complet, comme dans Empaqueter les artefacts :

import joblib

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

Le résultat :

 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

Ce que cela nous apprend. Le journal contient neuf lignes et la raison de chaque décision. Avec le pipeline enregistré, il permet à un autre analyste de reproduire le travail. Il permet aussi de réentraîner le modèle le mois prochain de la même façon.

Recommandations

Daniel écrit une note d’une page pour Amara, dans ses mots à elle, pas dans le langage des modèles :

Utilisez la prévision pour planifier les heures achalandées des jours de travail. Elle se trompe d’environ 45 vélos par heure en automne. C’est la moitié de l’erreur du profil sur 8 semaines que vous utilisez aujourd’hui.
Planifiez avec la prévision × 1,26 pour toute heure dont la prévision atteint 100 ou plus. En été, cela couvrait 9 heures achalandées sur 10. En automne, cela en couvrait 8 sur 10. Au début de chaque mois, nous vérifierons donc le facteur sur les 4 dernières semaines.
Attendez-vous à une prévision trop basse, d’environ 15 vélos par heure, parce que le réseau continue de croître. Nous réentraînons le modèle chaque mois avec les données les plus récentes.
Ajoutez les événements spéciaux à la main. Le modèle ne connaît pas les feux d’artifice, les courses ou les marches. Ces jours-là, votre expérience fait mieux que le modèle.
Les jours d’alertes météorologiques graves, décidez vous-même. Le modèle connaît la pluie, pas un ouragan ni une ville à l’arrêt.

Voici aussi à quoi le modèle ne doit pas servir. Il ne vous dit pas quelle station manquera de vélos, car les données couvrent seulement le réseau entier. Il ne sert pas à planifier combien de vélos acheter dans les prochaines années. Il n’a vu que deux années de croissance rapide. Il ne prévoit pas plus d’un jour à l’avance. Il ne doit servir à aucune décision concernant un membre du personnel en particulier.

Limites et prochaines étapes

Chaque résultat ci-dessus a des limites. Notez-les avant que quelqu’un d’autre les découvre :

Prévisions météo, pas météo observée. Les données contiennent la météo qui s’est réellement produite. Avec une vraie prévision faite la veille au soir, les erreurs seront plus grandes, surtout les jours de pluie.
Le réseau entier, pas les stations. Le rééquilibrage nécessite une prévision par station ou par secteur. Il faut les données des trajets de chaque station, que ce jeu de données ne contient pas.
Deux années de croissance rapide. Le modèle apprend un niveau qui change constamment. Le réentraînement mensuel et la vérification mensuelle du biais et du facteur de planification font partie du modèle. Ce ne sont pas des ajouts facultatifs.
Heures manquantes. Il manque les lignes de 165 heures, surtout pendant des fermetures. Le modèle n’a jamais appris à reconnaître une fermeture, donc il ne peut pas en prévoir une.
Valeurs inhabituelles des capteurs. La vitesse du vent vaut 0 pendant 2 180 heures, et l’humidité pendant 22 heures. Une prochaine étape consiste à les remplacer par des valeurs d’une autre source météo.

Voici la fiche du modèle, en bref, comme dans Documenter les limites d’utilisation :

Section	Contenu
Utilisation prévue	Prévisions à 20 h des locations horaires du réseau entier pour le lendemain, pour planifier le personnel et les fourgonnettes.
Utilisateurs	L’équipe des opérations.
Utilisations exclues	Rééquilibrage par station, planification à long terme, décisions sur un membre du personnel, météo extrême ou arrêt des activités de la ville.
Données	UCI Bike Sharing, Washington, D.C., 2011 et 2012, CC BY 4.0. Entraînement jusqu’en juin 2012, validation de juillet à septembre 2012, test d’octobre à décembre 2012.
Modèle	Amplification par gradient avec perte de Poisson, taux d’apprentissage de 0,05, 400 arbres, 11 variables calendaires et météorologiques.
Évaluation	MAE de validation de 38,7. Test (une seule fois) : MAE de 44,5, RMSE de 69,6, sous-estimation de 15,3. Référence de base du profil horaire sur le test : MAE de 94,1.
Sous-groupes	Heures de pointe : MAE de test de 88,7, sous-estimation de 39,9. Pluie légère : environ deux fois l’erreur par temps clair sur les données de validation. Sandy : MAE de 79,8.
Point de fonctionnement	Facteur de planification de 1,26 pour les prévisions de 100 ou plus; couverture de 0,891 en validation, de 0,788 au test.
Limites	Météo observée; réseau entier seulement; croissance; jours fériés avec événements.
Responsable et révision	Daniel Kim; révision mensuelle du biais, de la MAE des 4 dernières semaines et de la couverture du facteur de planification. Réentraîner quand la couverture des 4 dernières semaines descend sous 0,85.

Faites-en un projet de portfolio

Cette étude de cas constitue un bon premier projet pour un portfolio public. Elle présente des données réelles sous licence ouverte, une vraie décision et un résultat honnête. Pour vous l’approprier, changez au moins une chose. Par exemple, choisissez une autre date de séparation, une question par station ou les données ouvertes d’une autre ville. Publiez ensuite un petit dépôt.

Quoi mettre dans le dépôt :

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

Le README est la page que la plupart des gens lisent. Un modèle de README :

# 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.

Comment le décrire en entrevue, en deux phrases : « J’ai prévu la demande horaire de vélopartage pour planifier le personnel du lendemain. Une séparation aléatoire réduisait l’erreur apparente de 30 %, alors j’ai testé selon le temps; le modèle a réduit de moitié l’erreur de planification, et le test a montré qu’une marge de sécurité estivale ne suffit pas en automne. »

Pour aller plus loin :

Utilisez une vraie prévision météo plutôt que la météo observée, par exemple à partir d’archives météo publiques. Mesurez l’augmentation de l’erreur.
Faites des prévisions par station avec les données publiques des trajets d’un réseau de vélopartage. Planifiez les fourgonnettes par secteur.
Remplacez le facteur de planification par un modèle qui prévoit directement le 90e centile : HistGradientBoostingRegressor(loss="quantile", quantile=0.9).
Ajoutez un calendrier des événements publics et vérifiez l’erreur ces jours-là.
Effectuez un test rétrospectif mensuel sur 2012 et tracez la MAE et le biais par mois.

Vérification des connaissances

Vérifiez votre compréhension

La séparation aléatoire de Lena donne une MAE de 27 vélos par heure, et la séparation temporelle de Daniel donne 38,7. Les deux utilisent le même modèle. Expliquez la différence et dites à quel nombre Amara devrait s'attendre.

Afficher la réponse

Avec une séparation aléatoire, chaque heure de test a ses voisines dans les données d’entraînement : l’heure précédente et l’heure suivante. Les mêmes semaines de croissance et de météo s’y trouvent aussi. Le modèle apprend le niveau des journées mêmes sur lesquelles on le teste. Amara fait ses prévisions à 20 h pour demain, alors qu’aucune heure de demain n’est connue. La séparation temporelle reproduit cette situation. Il faut donc s’attendre à 38,7 en validation et à 44,5 sur le trimestre de test suivant. Même ces résultats sont un peu optimistes, car les données contiennent la météo observée, pas les prévisions.

Vérifiez votre compréhension

Amara demande l'erreur en pourcentage, car « une erreur de 45 vélos » ne dit rien au conseil d'administration. Quel problème pose la MAPE sur ces données, et que rapportez-vous à la place?

Afficher la réponse

La MAPE divise chaque erreur par le nombre réel. La nuit, une heure peut compter 1 ou 2 locations. Une erreur de 3 vélos devient alors 150 % ou 300 %. Ces heures dominent ensuite la moyenne. La MAPE indique donc une très mauvaise prévision, alors qu’elle est bonne là où cela compte. Rapportez plutôt la MAE comme part de la demande moyenne : 44,5 sur 219,2 locations par heure, soit environ 20 % au test. Ou rapportez la MAPE seulement pour les heures achalandées, en précisant laquelle vous avez utilisée.

Vérifiez votre compréhension

Sur le trimestre de test, le facteur de planification de 1,26 n'a couvert que 79 % des heures achalandées au lieu de 90 %. Daniel veut choisir un nouveau facteur sur ce trimestre. Que lui conseillez-vous?

Afficher la réponse

Ne le choisissez pas sur le trimestre de test. Cela transformerait le test en deuxième jeu de données de validation, et les 79 % ne seraient plus un résultat honnête. Le constat est que le facteur change avec la saison. Écrivez une règle. Au début de chaque mois, choisissez le facteur qui couvre 9 heures achalandées sur 10 dans les 4 dernières semaines de données réelles. À la fin du mois, vérifiez sa couverture. En utilisation réelle, les semaines les plus récentes jouent le rôle des données de validation.

Récapitulatif

Le problème a été cadré à partir de la décision : une heure du réseau entier, prévue à 20 h la veille. Les variables casual et registered sont donc des fuites de données.
La séparation suivait le temps, car le modèle prévoit l’avenir. Une séparation aléatoire donnait une MAE de 27, un résultat qu’Amara n’obtiendrait jamais.
Le profil horaire des 8 dernières semaines était une référence de base solide (MAE de validation de 52,2). Les modèles linéaires ne l’ont dépassée qu’après avoir appris que les nombres se multiplient et que chaque type de journée a sa forme.
L’amplification par gradient avec une perte de Poisson a gagné (MAE de validation de 38,7, jamais sous zéro). Les locations passées ont dégradé le résultat, et le journal conserve ce constat.
Un facteur de planification de 1,26 a transformé la prévision en plan. Sur le trimestre de test, utilisé une seule fois, le modèle se trompait de 44,5 vélos par heure, contre 94,1 pour le profil. Mais le facteur ne couvrait que 79 % des heures achalandées. Vérifiez donc les marges chaque mois.

Ensuite, une deuxième étude de cas de régression : la première estimation du prix d’une maison par un prêteur. La distribution des prix est asymétrique, les erreurs sont des pourcentages, et le marché évolue pendant que le modèle reste inchangé.

Leçon suivante : Estimer le prix des maisons