Fitting Poisson Match Models via MLE

Overview

The poisson-distribution page gives the probability machinery — this page gives the fitting recipe. The standard approach is Maher (1982): every team i gets an attack parameter αᵢ and a defense parameter βᵢ, and the expected goals for a match between home team i and away team j are

  • λ_home = αᵢ · βⱼ · γ (γ = home advantage multiplier)
  • λ_away = αⱼ · βᵢ

All parameters are estimated jointly over the full match history by maximizing the Poisson log-likelihood. This is fundamentally different from per-match fitting: a single optimization over ~2N+1 parameters (N teams) and thousands of matches, where every match informs both teams' parameters.

Dixon & Coles (1997) extended this with the low-score correction (see dixon-coles-correction) and — just as importantly — time-decay weighting, so older matches count less.

Why It Matters

This is the part a v1 implementation actually has to get right. Three details are easy to miss:

  1. Identifiability. The likelihood is invariant to scaling all attacks up and all defenses down (αᵢ → cαᵢ, βᵢ → βᵢ/c gives identical λs). Without a constraint the optimizer wanders along a ridge. Standard fix: constrain the mean log-attack to zero, Σ log αᵢ = 0 (or fix one team's α to 1).
  2. Joint estimation, not per-team averages. Attack/defense multipliers computed from raw goals-for/goals-against averages ignore schedule strength. MLE automatically adjusts for opponent quality because each match's likelihood contains the opponent's parameter.
  3. Time decay. Team strength drifts. Weighting match m's log-likelihood contribution by w_m = exp(−ξ · Δt_m) (Δt in days) keeps estimates current. ξ is a hyperparameter — tune it on out-of-sample predictive log-loss via walk-forward-validation, not on the training likelihood (the likelihood is always maximized at ξ=0).

Key Formula

Weighted log-likelihood over matches m = 1..M with home team h(m), away team a(m), score (x_m, y_m):

$$\ell = \sum_m w_m \Big[ \log P(x_m;\, \lambda_{h}) + \log P(y_m;\, \lambda_{a}) \Big], \quad w_m = e^{-\xi \, \Delta t_m}$$

with λ_h = α_{h(m)} β_{a(m)} γ and λ_a = α_{a(m)} β_{h(m)}, subject to Σᵢ log αᵢ = 0.

Optimize in log-space (work with log α, log β, log γ) so positivity is automatic and the constraint is linear. Typical ξ for international football: 0.001–0.005 per day (half-life of roughly 5 months to 2 years); club football tunes lower because teams play more often.

Code Snippet

import numpy as np
from scipy.optimize import minimize
from scipy.stats import poisson

def fit_poisson_mle(matches, teams, xi=0.002, ref_date=None):
    """Fit Maher-style Poisson model.

    matches: list of (date, home_team, away_team, home_goals, away_goals, neutral)
    teams:   list of team names
    Returns dict with per-team attack/defense and home advantage.
    """
    idx = {t: k for k, t in enumerate(teams)}
    n = len(teams)
    ref_date = ref_date or max(m[0] for m in matches)

    rows = []
    for date, ht, at, hg, ag, neutral in matches:
        w = np.exp(-xi * (ref_date - date).days)
        rows.append((idx[ht], idx[at], hg, ag, w, neutral))

    def neg_ll(params):
        log_a = params[:n] - params[:n].mean()   # enforce sum(log alpha) = 0
        log_b = params[n:2*n]
        log_gamma = params[-1]
        ll = 0.0
        for h, a, hg, ag, w, neutral in rows:
            adv = 0.0 if neutral else log_gamma   # no home advantage at neutral venues
            lam_h = np.exp(log_a[h] + log_b[a] + adv)
            lam_a = np.exp(log_a[a] + log_b[h])
            ll += w * (poisson.logpmf(hg, lam_h) + poisson.logpmf(ag, lam_a))
        return -ll

    x0 = np.zeros(2 * n + 1)
    res = minimize(neg_ll, x0, method="L-BFGS-B")
    log_a = res.x[:n] - res.x[:n].mean()
    return {
        "attack":  {t: np.exp(log_a[idx[t]]) for t in teams},
        "defense": {t: np.exp(res.x[n + idx[t]]) for t in teams},
        "home_advantage": np.exp(res.x[-1]),
        "converged": res.success,
    }

For speed on large datasets, vectorize the inner loop with numpy arrays indexed by team — the per-match Python loop is the bottleneck. Adding the Dixon-Coles ρ is one extra scalar parameter and a τ factor in the per-match likelihood.

Pitfalls

  • Forgetting the constraint — symptoms: optimizer reports success but parameters differ run-to-run while predictions are identical, or L-BFGS hits max iterations on a flat ridge.
  • Tuning ξ in-sample — the training likelihood always prefers ξ=0. Tune on held-out predictive log-loss only.
  • Teams with few matches — minnows that rarely appear get wild parameters. Shrink toward the mean (add a Gaussian prior on log α, log β — i.e. L2 penalty), or pool weak teams. Critical for international data; see international-football-modelling.
  • Neutral venues — apply γ only when there is a genuine home team (the neutral flag above). At the World Cup that means hosts only.
  • Separate home/away parameters — some implementations give each team distinct home and away attack/defense. That doubles parameters; with sparse international data, don't.

See Also