Poisson regression: when y is a count, a straight line predicts 1.4 calls — a log-link never can, and e^β is a rate ratio

The moment your target is a count — calls, defects, goals, claims, website visits — a straight line quietly lies about it. Counts are non-negative integers, they're skewed when small, and their variance grows with their mean. Fit ordinary least squares and it will happily predict −1.4 calls and assume the noise is a fixed symmetric bell — both wrong. Poisson regression is the generalized linear model built for exactly this, and I built a demo that fits it live in the browser with plain-JS gradient ascent on data you can place yourself. Here's what it taught me.

The model: a log-link that can't go negative

Assume each observation is drawn yᵢ ~ Poisson(λᵢ). The linear predictor η = β₀ + β₁x can be any real number — so you put a log-link on the mean, log λ = η, and the inverse λ = e^η squashes it into (0, ∞). The expected count is always a valid positive rate, no matter what the coefficients are.

def eta(beta, X):  return X @ beta               # β₀ + β₁x, any real number
def rate(beta, X): return np.exp(eta(beta, X))   # λ = e^η  (always > 0)

Fitting = maximizing the Poisson log-likelihood

There's no least-squares here. You maximize the probability the model assigns to the observed counts:

def log_likelihood(beta, X, y):
    lam = rate(beta, X)
    return np.sum(y * np.log(lam) - lam - gammaln(y + 1))

The log(y!) term is constant in β, so it doesn't move the optimum — it's there only so the number is a true log-likelihood.

The gradient collapses to one beautiful formula

Differentiate that and the messy chain rule cancels perfectly, thanks to the log-link, leaving the "score":

def gradient(beta, X, y):
    lam = rate(beta, X)
    return X.T @ (y - lam)          # ∇LL = Xᵀ(y − λ)

That's the sum of each feature row weighted by the residual — observed count minus predicted count. Setting it to zero is the maximum-likelihood condition, and it forces Σyᵢ = Σλᵢ: fitted totals match observed totals. The striking thing is that linear and logistic regression have the exact same gradient shape Xᵀ(y − μ) — pick the distribution of y, pick the link that maps β·x into the right range, and you've described the whole GLM family.

Climb it — the likelihood is concave, so there's one peak

Gradient ascent walks straight uphill to the single global maximum. Two touches make it bullet-proof on hand-drawn data: standardize x so the coordinates share a scale, and divide the step by the diagonal curvature Σλ·x² — because λ = e^η can explode, a raw fixed learning rate overflows, but this self-scaling keeps every step stable.

for _ in range(iters):
    lam  = np.exp(Z @ beta)         # λ = e^(Zβ)
    grad = Z.T @ (y - lam)          # score
    H    = (Z * Z).T @ lam          # diagonal curvature Σλ·z²
    beta += lr * grad / (H + 1e-9)  # curvature-scaled ascent

The payoff is interpretability

Because the model is linear on the log scale, coefficients act multiplicatively on the count. A one-unit rise in x multiplies the expected count by e^β₁ — the rate ratio. e^β₁ = 1.28 means "+1 in x ⇒ 28% more events." And when you're modelling a rate rather than a raw count — events per unit of exposure, like claims per policy-year — you add log(exposure) as a fixed offset column: log λ = log(exposure) + β·x. One caveat I flag in the demo: Poisson assumes variance equals the mean. Real counts are often overdispersed (variance ≫ mean) — check deviance ÷ degrees-of-freedom, and if it's well above 1, reach for the Negative Binomial (it adds a dispersion parameter); if zeros pile up, a zero-inflated model.

The takeaway I now carry: the moment y is a count, don't fit a line to it. Put a log-link on the mean and the geometry finally matches the data — the fitted curve can never dive below zero into the impossible band where OLS happily goes.

Drop some count data, hit Refit, and watch the coefficients climb the log-likelihood live:
https://dev48v.infy.uk/ml/day47-poisson-regression.html