← Retour à l'épreuve

Sweden SAIO · Correction IOAI Training

Relay the delay!

Cette page donne toute la méthode. Lisez-la une fois que vous avez cherché.

Métrique
MSE (minutes^2)
Évaluée sur
held-out validation (official test labels are private)
Notre score
val MSE 81.45 vs 87.15 const (-6.5 %)
Niveau attendu
MSE 80 for 100 points, 90 for 0

Notre solution n'atteint pas le niveau attendu : il reste de la marge pour faire mieux

Nous avons résolu l'épreuve nous-mêmes : téléchargé les données, mené les expériences et mesuré le score ci-dessous. L'énoncé et les données appartiennent à leurs organisateurs. Cette correction est rédigée en anglais.

The problem in brief

You get one million German train arrivals (July 2024 to January 2026, open data from Deutsche Bahn) with the station (anonymised as a number from 0 to 99), the scheduled arrival time and the delay in minutes. For 100 000 more arrivals from the same period you only know the station and the time, and you must predict the delay (an integer, possibly negative: trains can be early). There is no public/private split and only 8 submissions are allowed, so all model selection has to happen on your own validation.

How you are scored

The metric is the mean squared error S=1n∑i(xi−yi)2S = \frac{1}{n}\sum_i (x_i - y_i)^2, mapped to points with a base of 90 and a best of 80:

Points=100⋅max⁡(0,min⁡(1,90−S90−80)).\text{Points} = 100\cdot\max\Big(0,\min\Big(1,\frac{90 - S}{90-80}\Big)\Big).

Predicting the overall mean gives about 91 (the variance of the delays), so 90 is essentially "the constant", and 80 means explaining 11 % of the variance. We hold out a random 10 % of the training rows as validation, which mimics the test (random arrivals from the same period). On our validation split the constant gives 87.15, a bit lower than 91 by chance: a few extreme delays move the MSE of 100 000 rows by several units, so we compare relative gains rather than absolute numbers.

First look at the data

Left: histogram of delays on a log scale; most trains arrive 0 to 5 minutes late, with a long tail up to two hours and more. Middle: share of the squared error coming from trains whose delay is more than 5 to 200 minutes away from the mean; trains more than 60 minutes off are 0.45 % of the rows but half of the squared error. Right: the mean delay per day over the 19 months, fluctuating between 2 and 8 minutes, with some bad days.

  • The median delay is 1 minute, 90 % of trains are less than 10 minutes late, but the tail is long (1 % beyond 43 minutes, the maximum 521).
  • The metric is dominated by rare disasters. Trains more than an hour off the mean are 0.45 % of the rows but make 49 % of the squared error. A model can only reduce the error where the big delays are somewhat predictable.
  • There are a few absurd values around −1439-1439 minutes: an arrival recorded just after midnight for a train scheduled just before, a clock wrap-around (24 h = 1440 min). They are noise for any model.
  • Some days are much worse than others for the whole network (right panel), and stations differ in their average delay. Test rows are random arrivals from the same 19 months, so the test contains the same bad days as the training set: the date is a legitimate and very informative feature here, because it identifies the day's weather, strike or incident.

A first baseline

The organisers' notebook predicts 0 everywhere (worse than the mean, since the mean delay is 3.9 minutes). Predicting the training mean gives 87.15 on our validation split. Simple group means, smoothed toward the global mean, show which information matters (validation MSE):

group mean byMSE
nothing (constant)87.15
weekday86.99
hour of the day86.59
calendar day86.47
station and calendar day85.33
station84.47
station and hour of the day83.14

Station identity and its daily rhythm are the strongest signals; the calendar day adds information too, but each of these explains only a few percent.

The key idea

The delay is a sum of effects that interact: a given station is bad at 17:00 on weekdays, a storm on one day hurts some regions more, a strike hits every station. Group means can capture one combination at a time; gradient-boosted trees learn all the useful combinations from the raw features at once. A tree splits on, for example, "station in this set, and minute of the day between 16:30 and 18:45, and day index in this range", and boosting adds hundreds of such trees, each one correcting the residuals of the previous ones.

Two details matter because the metric is MSE.

  • Predict the mean, not the median. With a squared loss the best constant per situation is the conditional mean. Delays are skewed, so the mean (3.9) sits well above the median (1): an "intuitive" median-based prediction would lose points.
  • Tame the outliers in training, not in the predictions. We clip the training targets to [−30,150][-30, 150] minutes, so the clock wrap-arounds and the few extreme values cannot distort the fit. This barely changes the conditional means elsewhere (MSE 81.37 with clipping versus 81.45 without, within noise), and it protects the model against fitting noise.

The full solution

  1. Features: station (as a categorical variable), hour, minute of the day, weekday and day index since 1 July 2024.
  2. LightGBM, squared loss: 600 trees, learning rate 0.03, 63 leaves, at least 200 trains per leaf (a strong regulariser: with so much noise, a leaf must average many trains), row and column subsampling 0.8, targets clipped to [−30,150][-30, 150].
  3. Fit on all training rows, predict the test rows, round to integers (rounding adds at most 1/12≈0.081/12 \approx 0.08 to the MSE on average).
m = lgb.LGBMRegressor(n_estimators=600, learning_rate=0.03, num_leaves=63,
                      min_child_samples=200, subsample=0.8, subsample_freq=1,
                      colsample_bytree=0.8)
m.fit(features(train), np.clip(y, -30, 150), categorical_feature=["station"])

Results

Validation MSE on a random 10 % of the training rows (the test labels are private):

approachvalidation MSEgain vs constant
predict 0 (organisers' notebook)about 102worse
constant mean87.150 %
smoothed mean by station and hour of the day83.144.6 %
LightGBM on 5 raw features (final)81.456.5 %
LightGBM + 7 out-of-fold target encodings (station-day, day-hour...)81.476.5 %
LightGBM + local windows (mean delay at the same station within 15 min to 6 h, and network-wide within 10 to 30 min)81.706.3 %

If the same 6.5 % relative gain holds on the test set, whose constant error is about 90 to 91, our MSE would be about 84 to 85, which is roughly 55 to 60 points: an estimate, not a measured score.

The two failed rows are instructive. Hand-made encodings of "how late were the trains at this station that day" seemed promising, but the trees already found the same information from the raw day index and station, so the engineered features added nothing. The local time windows were even slightly worse: with only one million trains spread over 100 stations and 580 days, a station sees about 17 arrivals a day in the data, and the delay of a train 30 minutes earlier says little about the next one. The remaining error is mostly the large, unpredictable delays: half of the squared error comes from 0.45 % of the trains, and nothing in "station + scheduled time" announces them.

What to remember

  • Find out where the squared error lives. When 0.45 % of the rows make half of the MSE, expect modest gains and focus on calibrating the mean, not on the typical case.
  • Under MSE, predict conditional means: skewed targets make the mean very different from the median.
  • Time can be a feature when train and test share the same period. A random split across the same days lets the model learn "bad days"; with a future test period this would be leakage-free but useless.
  • Gradient boosting learns interactions (station x time of day x day) that one-at-a-time group means cannot.
  • Check that feature engineering adds information, not just a new view of features the model already has; and with only 8 submissions, trust a clean validation split.
Code complet (solution.py)
"""SAIO "Relay the delay!": predict train delays (minutes) from station and scheduled time.

Model: gradient-boosted trees (LightGBM, squared loss) on five raw features:
    station (categorical), hour, minute of the day, weekday, day index since 2024-07-01.
The day index lets the trees learn "bad days" (storms, strikes) shared by the whole
network, and station x minute-of-day the regular rush-hour patterns.
Training targets are clipped to [-30, 150] minutes: the -1439 values are clock
wrap-arounds, and a handful of extreme delays would otherwise dominate the fit.

Run: python solution.py --val     (90/10 random split of train.csv, prints MSE)
     python solution.py           (fits on all of train.csv, writes answer.txt)
"""
import sys
import numpy as np
import pandas as pd
import lightgbm as lgb

T0 = pd.Timestamp("2024-07-01")


def features(df):
    t = pd.to_datetime(df.time)
    return pd.DataFrame({
        "station": df.station_name.values,
        "h": t.dt.hour.values,
        "mod": (t.dt.hour * 60 + t.dt.minute).values,
        "wd": t.dt.weekday.values,
        "day": (t - T0).dt.days.values,
    })


def model():
    return lgb.LGBMRegressor(n_estimators=600, learning_rate=0.03, num_leaves=63,
                             min_child_samples=200, subsample=0.8, subsample_freq=1,
                             colsample_bytree=0.8, n_jobs=2, verbose=-1)


def score(m):
    return 100 * max(0.0, min(1.0, (90 - m) / 10))


if __name__ == "__main__":
    d = pd.read_csv("train.csv")
    y = d.delay.values.astype(float)
    if "--val" in sys.argv:
        hold = np.random.default_rng(0).random(len(d)) < 0.1
        Xa, Xb = features(d[~hold]), features(d[hold])
        m = model().fit(Xa, np.clip(y[~hold], -30, 150), categorical_feature=["station"])
        p = m.predict(Xb)
        const = ((y[hold] - y[~hold].mean()) ** 2).mean()
        mse = ((p - y[hold]) ** 2).mean()
        print(f"constant: {const:.2f}   GBM: {mse:.2f}   relative gain {100*(1-mse/const):.1f} %")
    else:
        t = pd.read_csv("test.csv")
        m = model().fit(features(d), np.clip(y, -30, 150), categorical_feature=["station"])
        p = np.rint(m.predict(features(t))).astype(int)    # the output must be integers
        open("answer.txt", "w").write("\n".join(map(str, p)) + "\n")