Skip to content

Bike-sharing forecasting

In this tutorial we're going to forecast the number of bikes in 5 bike stations from the city of Toulouse. We'll do so by building a simple model step by step. The dataset contains 182,470 observations. Let's first take a peak at the data.

from pprint import pprint
from river import datasets

X_y = datasets.Bikes()

for x, y in X_y:
    pprint(x)
    print(f'Number of available bikes: {y}')
    break

{'clouds': 75,
 'description': 'light rain',
 'humidity': 81,
 'moment': datetime.datetime(2016, 4, 1, 0, 0, 7),
 'pressure': 1017.0,
 'station': 'metro-canal-du-midi',
 'temperature': 6.54,
 'wind': 9.3}
Number of available bikes: 1

Let's start by using a simple linear regression on the numeric features. We can select the numeric features and discard the rest of the features using a Select. Linear regression is very likely to go haywire if we don't scale the data, so we'll use a StandardScaler to do just that. We'll evaluate the model by measuring the mean absolute error. Finally we'll print the score every 20,000 observations. 

from river import compose
from river import linear_model
from river import metrics
from river import evaluate
from river import preprocessing
from river import optim

X_y = datasets.Bikes()

model = compose.Select('clouds', 'humidity', 'pressure', 'temperature', 'wind')
model |= preprocessing.StandardScaler()
model |= linear_model.LinearRegression(optimizer=optim.SGD(0.001))

metric = metrics.MAE()

evaluate.progressive_val_score(X_y, model, metric, print_every=20_000)

[20,000] MAE: 4.912727
[40,000] MAE: 5.333554
[60,000] MAE: 5.330948
[80,000] MAE: 5.392313
[100,000] MAE: 5.423059
[120,000] MAE: 5.541223
[140,000] MAE: 5.613023
[160,000] MAE: 5.622428
[180,000] MAE: 5.567824





MAE: 5.563893

The model doesn't seem to be doing that well, but then again we didn't provide a lot of features. Generally, a good idea for this kind of problem is to look at an average of the previous values. For example, for each station we can look at the average number of bikes per hour. To do so we first have to extract the hour from the moment field. We can then use a TargetAgg to aggregate the values of the target.

from river import feature_extraction
from river import stats

X_y = iter(datasets.Bikes())

def get_hour(x):
    x['hour'] = x['moment'].hour
    return x

model = compose.Select('clouds', 'humidity', 'pressure', 'temperature', 'wind')
model += (
    get_hour |
    feature_extraction.TargetAgg(by=['station', 'hour'], how=stats.Mean())
)
model |= preprocessing.StandardScaler()
model |= linear_model.LinearRegression(optimizer=optim.SGD(0.001))

metric = metrics.MAE()

evaluate.progressive_val_score(X_y, model, metric, print_every=20_000)

[20,000] MAE: 3.721246
[40,000] MAE: 3.829972
[60,000] MAE: 3.845068
[80,000] MAE: 3.910259
[100,000] MAE: 3.888652
[120,000] MAE: 3.923727
[140,000] MAE: 3.980953
[160,000] MAE: 3.950034
[180,000] MAE: 3.934545





MAE: 3.933498

By adding a single feature, we've managed to significantly reduce the mean absolute error. At this point you might think that the model is getting slightly complex, and is difficult to understand and test. Pipelines have the advantage of being terse, but they aren't always to debug. Thankfully river has some ways to relieve the pain.

The first thing we can do it to visualize the pipeline, to get an idea of how the data flows through it.

model


['clouds', 'humidity', 'pressure', 'temperature', 'wind']
 {'whitelist': {'pressure', 'temperature', 'humidity', 'wind', 'clouds'}}

get_hour
 def get_hour(x): x['hour'] = x['moment'].hour return x

target_mean_by_station_and_hour
 {'by': ['station', 'hour'], 'feature_name': 'target_mean_by_station_and_hour', 'groups': defaultdict(functools.partial(<function deepcopy at 0x7ffeb9fb5700>, Mean: 0.), {'metro-canal-du-midi_0': Mean: 7.93981, 'metro-canal-du-midi_1': Mean: 8.179704, 'metro-canal-du-midi_10': Mean: 12.486815, 'metro-canal-du-midi_11': Mean: 11.675479, 'metro-canal-du-midi_12': Mean: 10.197409, 'metro-canal-du-midi_13': Mean: 10.650855, 'metro-canal-du-midi_14': Mean: 11.109123, 'metro-canal-du-midi_15': Mean: 11.068934, 'metro-canal-du-midi_16': Mean: 11.274958, 'metro-canal-du-midi_17': Mean: 8.459136, 'metro-canal-du-midi_18': Mean: 7.587469, 'metro-canal-du-midi_19': Mean: 7.734677, 'metro-canal-du-midi_2': Mean: 8.35824, 'metro-canal-du-midi_20': Mean: 7.582465, 'metro-canal-du-midi_21': Mean: 7.190665, 'metro-canal-du-midi_22': Mean: 7.486895, 'metro-canal-du-midi_23': Mean: 7.840791, 'metro-canal-du-midi_3': Mean: 8.656051, 'metro-canal-du-midi_4': Mean: 8.868445, 'metro-canal-du-midi_5': Mean: 8.99656, 'metro-canal-du-midi_6': Mean: 9.09966, 'metro-canal-du-midi_7': Mean: 8.852642, 'metro-canal-du-midi_8': Mean: 12.66712, 'metro-canal-du-midi_9': Mean: 13.412186, 'place-des-carmes_0': Mean: 4.720696, 'place-des-carmes_1': Mean: 3.390295, 'place-des-carmes_10': Mean: 8.575303, 'place-des-carmes_11': Mean: 9.321546, 'place-des-carmes_12': Mean: 10.511931, 'place-des-carmes_13': Mean: 11.392745, 'place-des-carmes_14': Mean: 10.735003, 'place-des-carmes_15': Mean: 10.198787, 'place-des-carmes_16': Mean: 9.941479, 'place-des-carmes_17': Mean: 9.125579, 'place-des-carmes_18': Mean: 7.660775, 'place-des-carmes_19': Mean: 6.847649, 'place-des-carmes_2': Mean: 2.232181, 'place-des-carmes_20': Mean: 9.626876, 'place-des-carmes_21': Mean: 11.602929, 'place-des-carmes_22': Mean: 10.405537, 'place-des-carmes_23': Mean: 7.700904, 'place-des-carmes_3': Mean: 1.371981, 'place-des-carmes_4': Mean: 1.051665, 'place-des-carmes_5': Mean: 0.984993, 'place-des-carmes_6': Mean: 2.039947, 'place-des-carmes_7': Mean: 3.850369, 'place-des-carmes_8': Mean: 3.792624, 'place-des-carmes_9': Mean: 5.957182, 'place-esquirol_0': Mean: 7.415789, 'place-esquirol_1': Mean: 5.244396, 'place-esquirol_10': Mean: 19.465005, 'place-esquirol_11': Mean: 22.976512, 'place-esquirol_12': Mean: 25.324159, 'place-esquirol_13': Mean: 25.428847, 'place-esquirol_14': Mean: 24.57762, 'place-esquirol_15': Mean: 24.416851, 'place-esquirol_16': Mean: 23.555125, 'place-esquirol_17': Mean: 22.062564, 'place-esquirol_18': Mean: 18.10623, 'place-esquirol_19': Mean: 11.916638, 'place-esquirol_2': Mean: 2.858635, 'place-esquirol_20': Mean: 13.346362, 'place-esquirol_21': Mean: 16.743318, 'place-esquirol_22': Mean: 15.562088, 'place-esquirol_23': Mean: 10.911134, 'place-esquirol_3': Mean: 1.155929, 'place-esquirol_4': Mean: 0.73306, 'place-esquirol_5': Mean: 0.668546, 'place-esquirol_6': Mean: 1.21265, 'place-esquirol_7': Mean: 3.107535, 'place-esquirol_8': Mean: 8.518696, 'place-esquirol_9': Mean: 15.470588, 'place-jeanne-darc_0': Mean: 6.541667, 'place-jeanne-darc_1': Mean: 5.99892, 'place-jeanne-darc_10': Mean: 5.735553, 'place-jeanne-darc_11': Mean: 5.616142, 'place-jeanne-darc_12': Mean: 5.787478, 'place-jeanne-darc_13': Mean: 5.817699, 'place-jeanne-darc_14': Mean: 5.657546, 'place-jeanne-darc_15': Mean: 6.224604, 'place-jeanne-darc_16': Mean: 5.796141, 'place-jeanne-darc_17': Mean: 5.743089, 'place-jeanne-darc_18': Mean: 5.674784, 'place-jeanne-darc_19': Mean: 5.833068, 'place-jeanne-darc_2': Mean: 5.598169, 'place-jeanne-darc_20': Mean: 6.015755, 'place-jeanne-darc_21': Mean: 6.242541, 'place-jeanne-darc_22': Mean: 6.141509, 'place-jeanne-darc_23': Mean: 6.493028, 'place-jeanne-darc_3': Mean: 5.180556, 'place-jeanne-darc_4': Mean: 4.779626, 'place-jeanne-darc_5': Mean: 4.67063, 'place-jeanne-darc_6': Mean: 4.611995, 'place-jeanne-darc_7': Mean: 4.960718, 'place-jeanne-darc_8': Mean: 5.552273, 'place-jeanne-darc_9': Mean: 6.249573, 'pomme_0': Mean: 3.301532, 'pomme_1': Mean: 2.312914, 'pomme_10': Mean: 11.791436, 'pomme_11': Mean: 12.976854, 'pomme_12': Mean: 13.962654, 'pomme_13': Mean: 11.692257, 'pomme_14': Mean: 11.180851, 'pomme_15': Mean: 11.939586, 'pomme_16': Mean: 12.267051, 'pomme_17': Mean: 12.132993, 'pomme_18': Mean: 11.399108, 'pomme_19': Mean: 6.37021, 'pomme_2': Mean: 2.144453, 'pomme_20': Mean: 5.279234, 'pomme_21': Mean: 6.254257, 'pomme_22': Mean: 6.568678, 'pomme_23': Mean: 5.235756, 'pomme_3': Mean: 1.563622, 'pomme_4': Mean: 0.947328, 'pomme_5': Mean: 0.924175, 'pomme_6': Mean: 1.287805, 'pomme_7': Mean: 1.299456, 'pomme_8': Mean: 2.94988, 'pomme_9': Mean: 7.89396}), 'how': Mean: 0., 'target_name': 'target'}

StandardScaler
 {'counts': Counter({'target_mean_by_station_and_hour': 182470, 'pressure': 182470, 'temperature': 182470, 'humidity': 182470, 'wind': 182470, 'clouds': 182470}), 'means': defaultdict(<class 'float'>, {'clouds': 30.315131254453505, 'humidity': 62.24244533347998, 'pressure': 1017.0563060996391, 'target_mean_by_station_and_hour': 9.468200635816528, 'temperature': 20.50980692716619, 'wind': 3.4184331122924543}), 'vars': defaultdict(<class 'float'>, {'clouds': 1389.0025610928221, 'humidity': 349.59967918503554, 'pressure': 33.298307526514115, 'target_mean_by_station_and_hour': 33.720872727055365, 'temperature': 34.70701720774977, 'wind': 4.473627075744674}), 'with_std': True}

LinearRegression
 {'_weights': {'target_mean_by_station_and_hour': 3.871921823842874, 'pressure': 2.148194846809029, 'temperature': -2.795077764828726, 'humidity': 3.8817929874945247, 'wind': -0.2612193604527132, 'clouds': -0.610647061945034}, '_y_name': None, 'clip_gradient': 1000000000000.0, 'initializer': Zeros (), 'intercept': 6.096564954881429, 'intercept_init': 0.0, 'intercept_lr': Constant({'learning_rate': 0.01}), 'l2': 0.0, 'loss': Squared({}), 'optimizer': SGD({'lr': Constant({'learning_rate': 0.001}), 'n_iterations': 182470})}

We can also use the debug_one method to see what happens to one particular instance. Let's train the model on the first 10,000 observations and then call debug_one on the next one. To do this, we will turn the Bike object into a Python generator with iter() function. The Pythonic way to read the first 10,000 elements of a generator is to use itertools.islice.

import itertools

X_y = iter(datasets.Bikes())

model = compose.Select('clouds', 'humidity', 'pressure', 'temperature', 'wind')
model += (
    get_hour |
    feature_extraction.TargetAgg(by=['station', 'hour'], how=stats.Mean())
)
model |= preprocessing.StandardScaler()
model |= linear_model.LinearRegression()

for x, y in itertools.islice(X_y, 10000):
    y_pred = model.predict_one(x)
    model.learn_one(x, y)

x, y = next(X_y)
print(model.debug_one(x))

0. Input
--------
clouds: 0 (int)
description: clear sky (str)
humidity: 52 (int)
moment: 2016-04-10 19:03:27 (datetime)
pressure: 1,001.00000 (float)
station: place-esquirol (str)
temperature: 19.00000 (float)
wind: 7.70000 (float)

1. Transformer union
--------------------
    1.0 Select
    ----------
    clouds: 0 (int)
    humidity: 52 (int)
    pressure: 1,001.00000 (float)
    temperature: 19.00000 (float)
    wind: 7.70000 (float)

    1.1 get_hour | target_mean_by_station_and_hour
    ----------------------------------------------
    target_mean_by_station_and_hour: 7.97175 (float)

clouds: 0 (int)
humidity: 52 (int)
pressure: 1,001.00000 (float)
target_mean_by_station_and_hour: 7.97175 (float)
temperature: 19.00000 (float)
wind: 7.70000 (float)

2. StandardScaler
-----------------
clouds: -1.36138 (float)
humidity: -1.73083 (float)
pressure: -1.26076 (float)
target_mean_by_station_and_hour: 0.05496 (float)
temperature: 1.76232 (float)
wind: 1.45841 (float)

3. LinearRegression
-------------------
Name                              Value      Weight     Contribution  
                      Intercept    1.00000    6.58252        6.58252  
                    temperature    1.76232    2.47030        4.35345  
                         clouds   -1.36138   -1.92255        2.61732  
target_mean_by_station_and_hour    0.05496    0.54167        0.02977  
                           wind    1.45841   -0.77720       -1.13348  
                       humidity   -1.73083    1.44921       -2.50833  
                       pressure   -1.26076    3.78529       -4.77234

Prediction: 5.16889

The debug_one method shows what happens to an input set of features, step by step.

And now comes the catch. Up until now we've been using the progressive_val_score method from the evaluate module. What this does it that it sequentially predicts the output of an observation and updates the model immediately afterwards. This way of doing is often used for evaluating online learning models, but in some cases it is the wrong approach.

The following paragraph is extremely important. When evaluating a machine learning model, the goal is to simulate production conditions in order to get a trust-worthy assessment of the performance of the model. In our case, we typically want to forecast the number of bikes available in a station, say, 30 minutes ahead. Then, once the 30 minutes have passed, the true number of available bikes will be available and we will be able to update the model using the features available 30 minutes ago. If you think about, this is exactly how a real-time machine learning system should work. The problem is that this isn't what the progressive_val_score method is emulating, indeed it is simply asking the model to predict the next observation, which is only a few minutes ahead, and then updates the model immediately. We can prove that this is flawed by adding a feature that measures a running average of the very recent values.

X_y = datasets.Bikes()

model = compose.Select('clouds', 'humidity', 'pressure', 'temperature', 'wind')
model += (
    get_hour |
    feature_extraction.TargetAgg(by=['station', 'hour'], how=stats.Mean()) + 
    feature_extraction.TargetAgg(by='station', how=stats.EWMean(0.5))
)
model |= preprocessing.StandardScaler()
model |= linear_model.LinearRegression()

metric = metrics.MAE()

evaluate.progressive_val_score(X_y, model, metric, print_every=20_000)

[20,000] MAE: 20.159286
[40,000] MAE: 10.458898
[60,000] MAE: 7.2759
[80,000] MAE: 5.715397
[100,000] MAE: 4.775094
[120,000] MAE: 4.138421
[140,000] MAE: 3.682591
[160,000] MAE: 3.35015
[180,000] MAE: 3.091398





MAE: 3.06414

The score we got is too good to be true. This is simply because the problem is too easy. What we really want is to evaluate the model by forecasting 30 minutes ahead and only updating the model once the true values are available. This can be done using the moment and delay parameters in the progressive_val_score method. The idea is that each observation of the stream of the data is shown twice to the model: once for making a prediction, and once for updating the model when the true value is revealed. The moment parameter determines which variable should be used as a timestamp, while the delay parameter controls the duration to wait before revealing the true values to the model.

import datetime as dt

model = compose.Select('clouds', 'humidity', 'pressure', 'temperature', 'wind')
model += (
    get_hour |
    feature_extraction.TargetAgg(by=['station', 'hour'], how=stats.Mean()) + 
    feature_extraction.TargetAgg(by='station', how=stats.EWMean(0.5))
)
model |= preprocessing.StandardScaler()
model |= linear_model.LinearRegression()

evaluate.progressive_val_score(
    dataset=datasets.Bikes(),
    model=model,
    metric=metrics.MAE(),
    moment='moment',
    delay=dt.timedelta(minutes=30),
    print_every=20_000
)

[20,000] MAE: 2.24812
[40,000] MAE: 2.240287
[60,000] MAE: 2.270287
[80,000] MAE: 2.28649
[100,000] MAE: 2.294264
[120,000] MAE: 2.275891
[140,000] MAE: 2.261411
[160,000] MAE: 2.285978
[180,000] MAE: 2.289353





MAE: 2.29304

The score we now have is much more realistic, as it is comparable with the related data science competition. Moreover, we can see that the model gets better with time, which feels better than the previous situations. The point is that progressive_val_score method can be used to simulate a production scenario, and is thus extremely valuable.

Now that we have a working pipeline in place, we can attempt to make it more accurate. As a simple example, we'll using a EWARegressor from the expert module to combine 3 linear regression model trained with different optimizers. The EWARegressor will run the 3 models in parallel and assign weights to each model based on their individual performance.

from river import expert
from river import optim

model = compose.Select('clouds', 'humidity', 'pressure', 'temperature', 'wind')
model += (
    get_hour |
    feature_extraction.TargetAgg(by=['station', 'hour'], how=stats.Mean())
)
model += feature_extraction.TargetAgg(by='station', how=stats.EWMean(0.5))
model |= preprocessing.StandardScaler()
model |= expert.EWARegressor([
    linear_model.LinearRegression(optim.SGD()),
    linear_model.LinearRegression(optim.RMSProp()),
    linear_model.LinearRegression(optim.Adam())
])

evaluate.progressive_val_score(
    dataset=datasets.Bikes(),
    model=model,
    metric=metrics.MAE(),
    moment='moment',
    delay=dt.timedelta(minutes=30),
    print_every=20_000
)

[20,000] MAE: 2.253263
[40,000] MAE: 2.242859
[60,000] MAE: 2.272001
[80,000] MAE: 2.287776
[100,000] MAE: 2.295292
[120,000] MAE: 2.276748
[140,000] MAE: 2.262146
[160,000] MAE: 2.286621
[180,000] MAE: 2.289925





MAE: 2.293604