Was 7–0 a freak? Checking a model against the scores it predicts

import csv
import random
from collections import Counter, defaultdict
from datetime import datetime
from math import exp, factorial, log

names = [f"{y % 100:02d}{(y + 1) % 100:02d}" for y in range(2000, 2026)]
HELD, TEST = names[16:21], names[21:]  # settings chosen on 2016/17 to 2020/21, checked once on 2021/22 to 2025/26

def season(s):  # (date, home, away, home goals, away goals, half-time home, half-time away) in date order
    with open(f"SC0_{s}.csv", encoding="latin-1") as f:
        games = [r for r in csv.DictReader(f) if r.get("FTR") in ("H", "D", "A")]
    games.sort(key=lambda r: datetime.strptime(r["Date"], "%d/%m/%Y" if len(r["Date"]) == 10 else "%d/%m/%y"))
    return [(r["Date"], r["HomeTeam"], r["AwayTeam"], int(r["FTHG"]), int(r["FTAG"]), int(r["HTHG"]), int(r["HTAG"]))
            for r in games]

data = {s: season(s) for s in names}

def at_least(lam, k):  # Poisson chance of k or more goals
    return 1 - sum(exp(-lam) * lam ** i / factorial(i) for i in range(k))

# Bayesian parts 4 and 5: before every match, each side's attack and defence from last season (worth k = 20 matches)
# updated with this season so far. One row per side per match: who, the model's expected goals, their spread, the score.
k, rows, rng = 20, [], random.Random(1)
for prev, s in zip(names, names[1:]):
    lf, la, ln = Counter(), Counter(), Counter()
    for _, h, a, x, y, *_ in data[prev]:
        lf[h] += x; la[h] += y; ln[h] += 1
        lf[a] += y; la[a] += x; ln[a] += 1
    mu = sum(x + y for *_, x, y, _, _ in data[prev]) / (2 * len(data[prev]))
    venue = {"home": sum(g[3] for g in data[prev]) / len(data[prev]), "away": sum(g[4] for g in data[prev]) / len(data[prev])}
    gf, ga, n = Counter(), Counter(), Counter()
    att = lambda t: (k * (lf[t] / ln[t] if ln[t] else 0.85 * mu) + gf[t], k + n[t])  # Gamma shape and rate
    dfn = lambda t: (k * (la[t] / ln[t] if ln[t] else 1.15 * mu) + ga[t], k + n[t])
    for date, h, a, x, y, hx, hy in data[s]:
        for side, them, where, goals in ((h, a, "home", x), (a, h, "away", y)):
            (sa, ra), (sd, rd) = att(side), dfn(them)
            lam = (sa / ra) * (sd / rd) / mu * venue[where] / mu
            # allowing for doubt: 100 plausible pairs of strengths from the Gamma posteriors
            doubt = [rng.gammavariate(sa, 1 / ra) * rng.gammavariate(sd, 1 / rd) / mu * venue[where] / mu for _ in range(100)]
            rows.append({"season": s, "date": date, "side": side, "them": them, "lam": lam, "doubt": doubt,
                         "venue": venue[where], "goals": goals})
        gf[h] += x; ga[h] += y; n[h] += 1
        gf[a] += y; ga[a] += x; n[a] += 1

print("Biggest wins since 2000/01:")
for date, h, a, x, y, *_ in sorted((g for s in names for g in data[s]), key=lambda g: -abs(g[3] - g[4]))[:5]:
    print(f"  {date} {h} {x}-{y} {a}")
print(f"  scores of 7 or more by one side: {sum(r['goals'] >= 7 for r in rows)}, in {len(rows)} side-matches, 2001/02 to 2025/26")

print("\nThe check, 2001/02 to 2025/26: real scores against what the model says should happen")
for label, test in (("no goals", lambda l: exp(-l)), ("5 or more", lambda l: at_least(l, 5)), ("7 or more", lambda l: at_least(l, 7))):
    real = sum((r["goals"] == 0) if label == "no goals" else (r["goals"] >= int(label[0])) for r in rows)
    plain = sum(test(r["lam"]) for r in rows)
    doubt = sum(sum(test(l) for l in r["doubt"]) / len(r["doubt"]) for r in rows)
    print(f"  {label}: real {real}, the model {plain:.0f}, the model allowing for doubt {doubt:.0f}")

print(f"\nGoals: real {sum(r['goals'] for r in rows)}, the model {sum(r['lam'] for r in rows):.0f}. By the model's expected goals:")
bands = defaultdict(list)
for r in rows:
    bands[min(int(r["lam"] / 0.5), 6) * 0.5].append(r)
for b, rs in sorted(bands.items()):
    print(f"  {b:.1f}{'+' if b == 3 else f' to {b + 0.5:.1f}'}: {len(rs)} side-matches, expected {sum(r['lam'] for r in rs) / len(rs):.2f} a game, "
          f"real {sum(r['goals'] for r in rs) / len(rs):.2f}")

ahead = defaultdict(list)  # do sides well ahead ease off? second-half goals by the leading side, by half-time lead
for s in names:
    for _, h, a, x, y, hx, hy in data[s]:
        for lead, later in ((hx - hy, x - hx), (hy - hx, y - hy)):
            if lead >= 0:
                ahead[min(lead, 4)].append(later)
print("Second-half goals by the side ahead (or level) at half-time:",
      ", ".join(f"{'4+' if d == 4 else d} up {sum(v) / len(v):.2f} ({len(v)})" for d, v in sorted(ahead.items())))

def shrunk(r, b):  # pull expected goals a share 1 - b of the way towards the average for that venue
    return r["venue"] + b * (r["lam"] - r["venue"])

def fit(rs, b):  # average Poisson log-likelihood: higher is better
    return sum(-shrunk(r, b) + r["goals"] * log(shrunk(r, b)) - log(factorial(r["goals"])) for r in rs) / len(rs)

held = [r for r in rows if r["season"] in HELD]
test = [r for r in rows if r["season"] in TEST]
b = max((i / 100 for i in range(40, 121)), key=lambda v: fit(held, v))
print(f"\nShrink chosen on the held-back seasons: keep {b:.2f} of each side's distance from the average")
for label, rs in (("held-back", held), ("test", test)):
    print(f"  {label} seasons, {len(rs)} side-matches, log-likelihood {fit(rs, 1):.4f} before, {fit(rs, b):.4f} after")
for label, check, count in (("no goals", lambda l: exp(-l), lambda g: g == 0), ("5 or more", lambda l: at_least(l, 5), lambda g: g >= 5),
                            ("7 or more", lambda l: at_least(l, 7), lambda g: g >= 7)):
    print(f"  test seasons, {label}: real {sum(count(r['goals']) for r in test)}, "
          f"model {sum(check(r['lam']) for r in test):.1f} before, {sum(check(shrunk(r, b)) for r in test):.1f} after")

print("\nThe 7-0s, and the chance the model gave the winners of scoring 7 or more:")
for r in rows:
    game = next(g for g in data[r["season"]] if g[0] == r["date"] and r["side"] in g[1:3])
    if r["goals"] == 7 and sorted(game[3:5]) == [0, 7]:
        print(f"  {r['date']} {game[1]} {game[3]}-{game[4]} {game[2]}: {r['side']} expected {shrunk(r, b):.2f}, "
              f"chance {at_least(shrunk(r, b), 7):.2%}, about 1 in {1 / at_least(shrunk(r, b), 7):.0f}")
expected = sum(at_least(shrunk(r, b), 7) for r in rows)
real = sum(r["goals"] >= 7 for r in rows)
low = sum(exp(-expected) * expected ** i / factorial(i) for i in range(real + 1))
print(f"Over 25 seasons the shrunk model expects {expected:.1f} scores of 7 or more, {expected / 25:.2f} a season; "
      f"there were {real}. Chance of {real} or fewer by luck: {low:.1%}")