Bike rental: a network

Examples

We continue with the data and the features of the previous page. The best Poisson regression needed an interaction between the hour and the day type that we found by looking at the data. A network can learn such interactions from the one-hot columns alone.

import numpy as np
import pandas as pd
import torch
import matplotlib.pyplot as plt
from sklearn.datasets import fetch_openml
from sklearn.compose import ColumnTransformer
from sklearn.preprocessing import OneHotEncoder, StandardScaler, FunctionTransformer
from sklearn.pipeline import Pipeline
from sklearn.model_selection import KFold, cross_val_score
from sklearn.metrics import mean_poisson_deviance

bikes = fetch_openml(data_id=42712, as_frame=True, parser="auto").frame
bikes = bikes[bikes["weather"] != "heavy_rain"]
y = bikes["count"].to_numpy(dtype=np.float32)
X = bikes.drop(columns=["count", "feel_temp"])

features = ColumnTransformer([
    ("onehot", OneHotEncoder(drop="first", sparse_output=False),
     ["hour", "weekday", "month", "year", "holiday", "workingday", "weather", "season"]),
    ("scale", StandardScaler(), ["temp", "humidity", "windspeed"])])
to_float32 = FunctionTransformer(lambda a: a.astype(np.float32))   # torch wants float32
print(features.fit_transform(X).shape)

(17376, 51)

The network

The output is the rate \(\lambda(x) = \exp(f_\theta(x))\) and the loss is the negative Poisson log-likelihood, as for the Poisson regression, with a network \(f_\theta\) instead of a linear function.

class PoissonNet(torch.nn.Module):
    def __init__(self, input_dim, hidden_dim=128, dropout=0.0):
        super().__init__()
        self.f = torch.nn.Sequential(
            torch.nn.Linear(input_dim, hidden_dim),
            torch.nn.Dropout(dropout),
            torch.nn.GELU(),
            torch.nn.BatchNorm1d(hidden_dim),
            torch.nn.Linear(hidden_dim, hidden_dim, bias=False),
            torch.nn.Dropout(dropout),
            torch.nn.GELU(),
            torch.nn.Linear(hidden_dim, 1))

    def forward(self, x):
        return torch.exp(self.f(x)).squeeze(-1)


class PoissonLoss(torch.nn.Module):
    def forward(self, rate, y):
        return torch.mean(rate - y * torch.log(rate + 1e-8))

1

A smooth variant of the relu, \(\mathrm{gelu}(a) = a\,\Phi(a)\) with \(\Phi\) the standard normal cumulative distribution function.

2

Batch normalization: every hidden neuron is standardized with the mean and the standard deviation of its activity in the current batch, then scaled and shifted by two learned parameters. In evaluation mode it uses running averages of these statistics instead.

3

Without bias, because a bias right after batch normalization would be redundant with its learned shift.

4

The rate, never negative.

5

The negative Poisson log-likelihood without the term \(\log y!\), which does not depend on the parameters. The small constant avoids \(\log 0\).

The learning curve

Before cross-validating, we train once on the first 80% of the hours and watch the loss on the last 20%. This is how we chose the width, the batch size, the learning rate and the number of epochs.

n_fit = int(0.8 * len(y))
pre = Pipeline([("features", features), ("float", to_float32)]).fit(X.iloc[:n_fit])
X_fit, X_val = torch.tensor(pre.transform(X.iloc[:n_fit])), torch.tensor(pre.transform(X.iloc[n_fit:]))
y_fit, y_val = torch.tensor(y[:n_fit]), torch.tensor(y[n_fit:])

torch.manual_seed(1)
net = PoissonNet(X_fit.shape[1])
opt = torch.optim.AdamW(net.parameters(), lr=1e-3)
loss_fn = PoissonLoss()
loader = torch.utils.data.DataLoader(
    torch.utils.data.TensorDataset(X_fit, y_fit), batch_size=64, shuffle=True)
curve = []
for epoch in range(30):
    net.train()
    for xb, yb in loader:
        opt.zero_grad()
        loss_fn(net(xb), yb).backward()
        opt.step()
    net.eval()
    with torch.no_grad():
        curve.append((mean_poisson_deviance(y_fit, net(X_fit)),
                      mean_poisson_deviance(y_val, net(X_val))))

fig, ax = plt.subplots(figsize=(6, 3.2))
ax.plot(range(1, 31), [c[0] for c in curve], label="training")
ax.plot(range(1, 31), [c[1] for c in curve], label="validation")
ax.set(xlabel="epoch", ylabel="Poisson deviance", ylim=(0, None))
ax.legend()
plt.show()

Figure 46.1: Poisson deviance on the first 80% of the hours, which we train on, and on the last 20%, after every epoch.

The training deviance keeps falling, the validation deviance levels off much higher.

Cross-validation with skorch

To cross-validate the network exactly like the regressions, with the preprocessing inside the folds, we wrap it with the skorch package. Its NeuralNetRegressor has the fit and predict methods of an sklearn model and runs the training loop above inside fit.

from skorch import NeuralNetRegressor


def poisson_net(**kwargs):
    net = NeuralNetRegressor(
        PoissonNet,
        module__input_dim=features.fit_transform(X).shape[1],
        criterion=PoissonLoss,
        optimizer=torch.optim.AdamW,
        lr=1e-3, batch_size=64, max_epochs=30,
        iterator_train__shuffle=True,
        train_split=None,
        verbose=0,
        **kwargs)
    return Pipeline([("features", features), ("float", to_float32), ("net", net)])


torch.manual_seed(1)
cv = KFold(5)
dev = -cross_val_score(poisson_net(), X, y, cv=cv, scoring="neg_mean_poisson_deviance")
print(f"Poisson network, CV deviance: {dev.mean():.1f}  (folds: {np.round(dev, 1)})")

1

Arguments that start with module__ go to PoissonNet.

2

By default skorch holds back a fifth of the training data for validation. We switch that off, because the cross-validation does the validation.

Poisson network, CV deviance: 18.9  (folds: [16.7 13.4 14.6 25.1 24.6])

On the previous page the regression with one-hot hours had about 38, and the one with an hour by day type interaction about 19. The network gets there from the one-hot columns without being told about the interaction.

torch.manual_seed(1)
full = poisson_net().fit(X, y)
print(f"training deviance: {mean_poisson_deviance(y, full.predict(X)):.1f}")

training deviance: 6.6

The training deviance is much lower than the cross-validated one, so the network overfits. Weight decay, the L2 regularization of week 7, is a hyper-parameter of the optimizer.

for wd in [1e-2, 1e-1]:
    torch.manual_seed(1)
    model = poisson_net(optimizer__weight_decay=wd)
    dev = -cross_val_score(model, X, y, cv=cv, scoring="neg_mean_poisson_deviance")
    print(f"weight decay {wd}: CV deviance {dev.mean():.1f}")

1

Arguments that start with optimizer__ go to AdamW.

weight decay 0.01: CV deviance 18.9
weight decay 0.1: CV deviance 16.0

Reading the residuals

from sklearn.model_selection import cross_val_predict
from sklearn.linear_model import PoissonRegressor

torch.manual_seed(1)
pred_net = cross_val_predict(poisson_net(), X, y, cv=cv)
glm = Pipeline([("features", features), ("glm", PoissonRegressor(alpha=1e-6, max_iter=5000))])
pred_glm = cross_val_predict(glm, X, y, cv=cv)

fig, axes = plt.subplots(1, 2, figsize=(9, 3.4))
axes[0].plot(y, pred_net, ".", ms=1.5, alpha=0.3)
axes[0].plot([0, 1000], [0, 1000], "k--", lw=1)
axes[0].set(xlabel="observed count", ylabel="predicted rate (network)")
for pred, label in [(pred_glm, "regression, one-hot hour"), (pred_net, "network")]:
    r = pd.Series(y - pred).groupby(X["hour"].to_numpy()).mean()
    axes[1].plot(r.index, r.values, "o-", ms=3, label=label)
axes[1].axhline(0, color="k", lw=1)
axes[1].set(xlabel="hour of day", ylabel="mean residual")
axes[1].legend()
plt.tight_layout()
plt.show()

Figure 46.2: Left, the cross-validated prediction of every hour against the observed count. Right, the mean residual (observed minus predicted) per hour of the day, for the network and for the regression with one-hot hours.

The residuals of the regression with one-hot hours vary systematically with the hour, which is the interaction it cannot express. Those of the network stay closer to zero.

The final model

The cross-validation chose the network with weight decay 0.1. For predictions on new hours, we fit it once more on all the data.

torch.manual_seed(1)
final = poisson_net(optimizer__weight_decay=1e-1).fit(X, y)
final.predict(X.iloc[:3])

array([22.04892, 18.40782, 12.58117], dtype=float32)