Skip to main content

A Poisson network

Examples

We keep the Poisson likelihood from week 3 and replace the linear predictor by a network,

\[ \lambda(x) = \exp\bigl(f_\theta(x)\bigr), \qquad \mathcal L(\theta) = \frac1n\sum_{i=1}^n \bigl(\lambda(x_i) - y_i\log\lambda(x_i)\bigr). \]

The second expression is the negative Poisson log-likelihood without the term that does not depend on \(\theta\).

The network

import numpy as np
import pandas as pd
import torch
import torch.nn as nn
import matplotlib.pyplot as plt
from sklearn.model_selection import train_test_split
from sklearn.preprocessing import StandardScaler

torch.manual_seed(8)
rng = np.random.default_rng(8)

n = 4000
hour = rng.integers(0, 24, n)
weekday = rng.integers(0, 7, n)
temp = rng.uniform(-4, 34, n)
humidity = rng.uniform(0.2, 1.0, n)
holiday = (rng.random(n) < 0.05).astype(float)

peak = np.where(weekday < 5,
                np.exp(-0.5 * ((hour - 8) / 1.3) ** 2) + np.exp(-0.5 * ((hour - 18) / 1.6) ** 2),
                0.9 * np.exp(-0.5 * ((hour - 14) / 4.0) ** 2))
rate = np.exp(2.0 + 1.6 * peak + 0.05 * temp - 1.3 * humidity - 0.4 * holiday)
count = rng.poisson(rate)

X = np.column_stack([np.sin(2 * np.pi * hour / 24), np.cos(2 * np.pi * hour / 24),
                     (weekday < 5).astype(float), temp, humidity, holiday])

Xtr, Xte, ytr, yte = train_test_split(X, count, test_size=0.25, random_state=0)
sc = StandardScaler().fit(Xtr)
Xtr_t = torch.tensor(sc.transform(Xtr), dtype=torch.float32)
Xte_t = torch.tensor(sc.transform(Xte), dtype=torch.float32)
ytr_t = torch.tensor(ytr, dtype=torch.float32)
yte_t = torch.tensor(yte, dtype=torch.float32)
1
The scaler is fitted on the training part only. 
class PoissonNet(nn.Module):
    def __init__(self, input_dim, hidden=64, dropout=0.1):
        super().__init__()
        self.fc = nn.Sequential(
            nn.Linear(input_dim, hidden),
            nn.BatchNorm1d(hidden),
            nn.GELU(),
            nn.Dropout(dropout),
            nn.Linear(hidden, hidden),
            nn.GELU(),
            nn.Dropout(dropout),
            nn.Linear(hidden, 1),
        )

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


def poisson_loss(rate, y):
    return torch.mean(rate - y * torch.log(rate + 1e-8))


model = PoissonNet(X.shape[1])
opt = torch.optim.AdamW(model.parameters(), lr=1e-2, weight_decay=1e-4)
1
The exponential is the inverse link. It keeps the rate positive.
2
The negative Poisson log-likelihood, up to a constant. The small term inside the log guards against a rate of exactly zero.
3
weight_decay is the L2 penalty.

Training

loader = torch.utils.data.DataLoader(
    torch.utils.data.TensorDataset(Xtr_t, ytr_t), batch_size=64, shuffle=True)

train_curve, test_curve = [], []
for epoch in range(40):
    model.train()
    for xb, yb in loader:
        opt.zero_grad()
        poisson_loss(model(xb), yb).backward()
        opt.step()
    model.eval()
    with torch.no_grad():
        train_curve.append(float(poisson_loss(model(Xtr_t), ytr_t)))
        test_curve.append(float(poisson_loss(model(Xte_t), yte_t)))

print(f"final training loss: {train_curve[-1]:.4f}")
print(f"final test loss:     {test_curve[-1]:.4f}")
1
Training mode. Dropout is active and batch norm uses the batch statistics.
2
Evaluation mode. Dropout is off and batch norm uses the running statistics. Forgetting this is a common bug, and it makes the model look worse than it is. 
final training loss: -39.2777
final test loss:     -38.7020
fig, ax = plt.subplots()
ax.plot(train_curve, label="training")
ax.plot(test_curve, label="test")
ax.set(xlabel="epoch", ylabel="Poisson loss")
ax.legend()
plt.show()
Figure 41.1: Training and test loss. If the test curve turned up we would stop earlier.

Comparing to the GLM

from sklearn.linear_model import PoissonRegressor
from sklearn.metrics import mean_poisson_deviance

glm = PoissonRegressor(alpha=1e-4, max_iter=5000).fit(sc.transform(Xtr), ytr)

model.eval()
with torch.no_grad():
    net_pred = model(Xte_t).numpy()
glm_pred = glm.predict(sc.transform(Xte))

print(f"Poisson GLM deviance:     {mean_poisson_deviance(yte, glm_pred):.4f}")
print(f"Poisson network deviance: {mean_poisson_deviance(yte, net_pred):.4f}")
Poisson GLM deviance:     5.4424
Poisson network deviance: 1.1480
fig, ax = plt.subplots()
ax.plot(yte, glm_pred, ".", alpha=0.35, label="Poisson GLM")
ax.plot(yte, net_pred, ".", alpha=0.35, label="Poisson network")
lim = [0, yte.max() * 1.05]
ax.plot(lim, lim, ls="--", lw=1, color="#7a838b")
ax.set(xlabel="true count", ylabel="predicted count", xlim=lim, ylim=lim)
ax.legend()
plt.show()
Figure 41.2: True counts against predictions on the test set.

If the network does not beat the GLM with good features, we keep the GLM. It is faster, it is easier to explain, and it has fewer ways to go wrong.

Reading the residuals

hour_te = np.round(np.arctan2(Xte[:, 0], Xte[:, 1]) * 24 / (2 * np.pi)) % 24
resid = yte - net_pred

fig, ax = plt.subplots()
means = pd.Series(resid).groupby(hour_te).mean()
ax.plot(means.index, means.values, "o-")
ax.axhline(0, color="#a8452f", lw=1.2, ls="--")
ax.set(xlabel="hour of day", ylabel="mean residual")
plt.show()
Figure 41.3: Residuals against the hour of the day. A pattern here would mean the hour features are not doing their job.

A pattern left in the residuals is a predictor we have not used yet, not randomness. This plot is worth making for every predictor.

Where it goes wrong

Four mistakes come back every year.

Standardizing before the split. The test mean leaks into the training set. Fit the scaler on the training part only.

Using squared error on counts. It fits the busy hours and ignores the quiet ones, and it predicts negative rentals at night.

Tuning against the test set. After the third look it is a validation set, and the number you report is optimistic.

Forgetting model.eval(). Dropout stays on and the predictions are noisy.