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
dataset = datasets.Bikes()
for x, y in dataset:
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
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(
dataset.take(50000), model, metric, print_every=5_000
)
[5,000] MAE: 4.258099 [10,000] MAE: 4.495612 [15,000] MAE: 4.752074 [20,000] MAE: 4.912727 [25,000] MAE: 4.934188 [30,000] MAE: 5.164331 [35,000] MAE: 5.320877 [40,000] MAE: 5.333554 [45,000] MAE: 5.354958 [50,000] MAE: 5.378699
MAE: 5.378699
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
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(1e-2))
metric = metrics.MAE()
evaluate.progressive_val_score(
dataset.take(50000), model, metric, print_every=5_000
)
[5,000] MAE: 69.042914 [10,000] MAE: 36.269638 [15,000] MAE: 25.241059 [20,000] MAE: 19.781737 [25,000] MAE: 16.605912 [30,000] MAE: 14.402878 [35,000] MAE: 12.857216 [40,000] MAE: 11.647737 [45,000] MAE: 10.646566 [50,000] MAE: 9.94726
MAE: 9.94726
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']
(
clouds
humidity
pressure
temperature
wind
)
get_hour
def get_hour(x):
x['hour'] = x['moment'].hour
return x
y_mean_by_station_and_hour
(
by=['station', 'hour']
how=Mean ()
target_name="y"
)
StandardScaler
(
with_std=True
)
LinearRegression
(
optimizer=SGD (
lr=Constant (
learning_rate=0.01
)
)
loss=Squared ()
l2=0.
l1=0.
intercept_init=0.
intercept_lr=Constant (
learning_rate=0.01
)
clip_gradient=1e+12
initializer=Zeros ()
)
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 is that it sequentially predicts the output of an observation and updates the model immediately afterwards. This way of proceeding is often used for evaluating online learning models. But in some cases it is the wrong approach.
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.
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 in 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
evaluate.progressive_val_score(
dataset=dataset.take(50000),
model=model.clone(),
metric=metrics.MAE(),
moment="moment",
delay=dt.timedelta(minutes=30),
print_every=5_000,
)
[5,000] MAE: 4.675207 [10,000] MAE: 4.352476 [15,000] MAE: 4.193511 [20,000] MAE: 4.203433 [25,000] MAE: 4.226929 [30,000] MAE: 4.191629 [35,000] MAE: 4.227425 [40,000] MAE: 4.195404 [45,000] MAE: 4.102599 [50,000] MAE: 4.117846
MAE: 4.117846
The performance is a bit worse, which is to be expected. Indeed, the task is more difficult: the model is only shown the ground truth 30 minutes after making a prediction.
from deep_river.regression import Regressor
from river import feature_extraction
from river import stats
import torch
class LinearRegression(torch.nn.Module):
def __init__(self, n_features, outputSize=1):
super(LinearRegression, self).__init__()
self.linear = torch.nn.Linear(n_features, outputSize)
def forward(self, x):
out = self.linear(x)
return out
model = compose.Select("clouds", "humidity", "pressure", "temperature", "wind")
model += get_hour | feature_extraction.TargetAgg(
by=["station", "hour"], how=stats.Mean()
)
model |= preprocessing.StandardScaler()
model |= Regressor(
module=LinearRegression,
loss_fn="mse",
optimizer_fn="sgd",
lr=1e-2,
)
import datetime as dt
evaluate.progressive_val_score(
dataset=dataset.take(50000),
model=model.clone(),
metric=metrics.MAE(),
moment="moment",
delay=dt.timedelta(minutes=30),
print_every=5_000,
)
[5,000] MAE: 4.686527 [10,000] MAE: 4.358232 [15,000] MAE: 4.197241 [20,000] MAE: 4.20616 [25,000] MAE: 4.229106 [30,000] MAE: 4.193436 [35,000] MAE: 4.228976 [40,000] MAE: 4.196765 [45,000] MAE: 4.103809 [50,000] MAE: 4.118935
MAE: 4.118935
Building RNN Models¶
from deep_river.regression import Regressor, RollingRegressor
from river import feature_extraction
from river import stats
import torch
class RnnModule(torch.nn.Module):
def __init__(self, n_features, hidden_size):
super().__init__()
self.n_features = n_features
self.rnn = torch.nn.RNN(
input_size=n_features, hidden_size=hidden_size, num_layers=1
)
self.fc = torch.nn.Linear(in_features=hidden_size, out_features=1)
def forward(self, X, **kwargs):
output, hn = self.rnn(X) # lstm with input, hidden, and internal state
return self.fc(output[-1, :])
model = compose.Select("clouds", "humidity", "pressure", "temperature", "wind")
model += get_hour | feature_extraction.TargetAgg(
by=["station", "hour"], how=stats.Mean()
)
model |= preprocessing.StandardScaler()
model |= RollingRegressor(
module=RnnModule,
loss_fn="mse",
optimizer_fn="sgd",
lr=1e-2,
hidden_size=20,
window_size=32,
)
import datetime as dt
evaluate.progressive_val_score(
dataset=dataset.take(50000),
model=model.clone(),
metric=metrics.MAE(),
moment="moment",
delay=dt.timedelta(minutes=30),
print_every=5_000,
)
[5,000] MAE: 4.546407 [10,000] MAE: 4.531478 [15,000] MAE: 4.596383 [20,000] MAE: 4.826214 [25,000] MAE: 4.83528 [30,000] MAE: 5.137052 [35,000] MAE: 5.341721 [40,000] MAE: 5.319061 [45,000] MAE: 5.238281 [50,000] MAE: 5.25946
MAE: 5.25946
Building LSTM Models¶
class LstmModule(torch.nn.Module):
def __init__(self, n_features, hidden_size=1):
super().__init__()
self.n_features = n_features
self.hidden_size = hidden_size
self.lstm = torch.nn.LSTM(
input_size=n_features,
hidden_size=hidden_size,
num_layers=1,
bidirectional=False,
)
self.fc = torch.nn.Linear(in_features=hidden_size, out_features=1)
def forward(self, X, **kwargs):
output, (hn, cn) = self.lstm(
X
) # lstm with input, hidden, and internal state
return self.fc(output[-1, :])
model = compose.Select("clouds", "humidity", "pressure", "temperature", "wind")
model += get_hour | feature_extraction.TargetAgg(
by=["station", "hour"], how=stats.Mean()
)
model |= preprocessing.StandardScaler()
model |= RollingRegressor(
module=LstmModule,
loss_fn="mse",
optimizer_fn="sgd",
lr=1e-2,
hidden_size=20,
window_size=32,
)
import datetime as dt
evaluate.progressive_val_score(
dataset=dataset.take(50000),
model=model.clone(),
metric=metrics.MAE(),
moment="moment",
delay=dt.timedelta(minutes=30),
print_every=5_000,
)
[5,000] MAE: 4.236382 [10,000] MAE: 3.981297 [15,000] MAE: 3.838141 [20,000] MAE: 3.838378 [25,000] MAE: 3.878212 [30,000] MAE: 3.827275 [35,000] MAE: 3.889832 [40,000] MAE: 3.848552 [45,000] MAE: 3.770227 [50,000] MAE: 3.79234
MAE: 3.79234