Clue by clue. Naive Bayes, and the danger of counting twice

import csv
from collections import Counter, defaultdict
from datetime import datetime
from math import exp, log, pi, sqrt
from statistics import correlation

POINTS = {"H": (3, 0), "D": (1, 1), "A": (0, 3)}
RESULTS = "HDA"
CLUES = ["recent form", "this season", "last season"]
names = [f"{y % 100:02d}{(y + 1) % 100:02d}" for y in range(2000, 2026)]

def season(s):
    with open(f"SC0_{s}.csv", encoding="latin-1") as f:
        games = [r for r in csv.DictReader(f) if r.get("FTR") in POINTS]
    games.sort(key=lambda r: datetime.strptime(r["Date"], "%d/%m/%Y" if len(r["Date"]) == 10 else "%d/%m/%y"))
    return games

def points_per_game(games):
    pts, n = Counter(), Counter()
    for r in games:
        for team, p in zip((r["HomeTeam"], r["AwayTeam"]), POINTS[r["FTR"]]):
            pts[team] += p
            n[team] += 1
    return {t: pts[t] / n[t] for t in n}

# the same three gaps as part 8, home side minus away side, in points a game, all known before kick-off
rows = []
for s_last, s in zip(names, names[1:]):
    last = points_per_game(season(s_last))
    promoted = 0.85 * sum(last.values()) / len(last)
    history = defaultdict(list)
    for r in season(s):
        h, a = r["HomeTeam"], r["AwayTeam"]
        if len(history[h]) >= 5 and len(history[a]) >= 5:
            rows.append({"season": s, "match": f"{h} {r['FTHG']}-{r['FTAG']} {a}", "date": r["Date"], "y": r["FTR"], "x": [
                (sum(history[h][-5:]) - sum(history[a][-5:])) / 5,
                sum(history[h]) / len(history[h]) - sum(history[a]) / len(history[a]),
                last.get(h, promoted) - last.get(a, promoted)]})
        for team, p in zip((h, a), POINTS[r["FTR"]]):
            history[team].append(p)

def with_strength(data, w):  # one combined clue: w of this season's gap plus (1 - w) of last season's
    for r in data:
        r["x"] = r["x"][:3] + [w * r["x"][1] + (1 - w) * r["x"][2]]

def naive_bayes(fit, clues):
    """Base rates, then for each clue a bell curve per result: its average and spread among home wins, draws, away wins."""
    prior = Counter(r["y"] for r in fit)
    curves = {}
    for c in RESULTS:
        xs = [r["x"] for r in fit if r["y"] == c]
        for j in clues:
            m = sum(x[j] for x in xs) / len(xs)
            curves[c, j] = (m, sqrt(sum((x[j] - m) ** 2 for x in xs) / len(xs)))
    def forecast(x, upto=None):  # upto: only the first few clues, to watch the odds move clue by clue
        score = {c: log(prior[c] / len(fit)) for c in RESULTS}
        for j in clues[:upto]:
            for c in RESULTS:
                m, sd = curves[c, j]
                score[c] += -log(sd * sqrt(2 * pi)) - (x[j] - m) ** 2 / (2 * sd ** 2)  # log of the bell curve's height
        top = max(score.values())
        e = {c: exp(v - top) for c, v in score.items()}
        return {c: e[c] / sum(e.values()) for c in RESULTS}
    return forecast, prior, curves

def logistic(fit):  # part 8's model, standardised and refitted here, for a match-by-match comparison
    mean = [sum(r["x"][j] for r in fit) / len(fit) for j in range(3)]
    sd = [sqrt(sum((r["x"][j] - mean[j]) ** 2 for r in fit) / len(fit)) for j in range(3)]
    z = lambda x: [(x[j] - mean[j]) / sd[j] for j in range(3)]
    def predict(w, v):
        e = {c: exp(w[c][0] + sum(a * b for a, b in zip(w[c][1:], v))) for c in RESULTS}
        return {c: e[c] / sum(e.values()) for c in RESULTS}
    w = {c: [0.0] * 4 for c in RESULTS}
    data = [(z(r["x"]), r["y"]) for r in fit]
    for _ in range(200):
        grad = {c: [0.0] * 4 for c in "HA"}
        for v, y in data:
            p = predict(w, v)
            for c in "HA":
                grad[c] = [g + (p[c] - (c == y)) * t for g, t in zip(grad[c], [1] + v)]
        for c in "HA":
            w[c] = [a - g / len(data) for a, g in zip(w[c], grad[c])]
    return lambda x: predict(w, z(x))

def log_loss(model, data):
    return sum(-log(model(r["x"])[r["y"]]) for r in data) / len(data)

def base_rates(fit):
    n = Counter(r["y"] for r in fit)
    return lambda x: {c: n[c] / len(fit) for c in RESULTS}

fit = [r for r in rows if r["season"] < "1617"]
held = [r for r in rows if "1617" <= r["season"] < "2122"]
train = [r for r in rows if r["season"] < "2122"]
test = [r for r in rows if r["season"] >= "2122"]
pairs = ((0, 1), (0, 2), (1, 2))
print("how alike the clues are (training matches): " + ", ".join(
    f"{CLUES[i]} and {CLUES[j]} {correlation([r['x'][i] for r in train], [r['x'][j] for r in train]):.2f}" for i, j in pairs))

# 1. on the held-back seasons: which clues, and how to combine the two season gaps
print(f"held back: base rates {log_loss(base_rates(fit), held):.4f}")
for clues in ([0, 1, 2], [0], [1], [2]):
    print(f"  naive Bayes on {', '.join(CLUES[j] for j in clues)}: {log_loss(naive_bayes(fit, clues)[0], held):.4f}")
scores = {}
for w in (0.3, 0.4, 0.5, 0.6, 0.7):
    with_strength(fit + held, w)
    scores[w] = log_loss(naive_bayes(fit, [3])[0], held)
    print(f"  one strength clue, {w:.0%} this season: {scores[w]:.4f}")
best = min(scores, key=scores.get)
print(f"chosen: {best:.0%} this season, {1 - best:.0%} last season")

# 2. refit on all training matches; score the test seasons once
with_strength(rows, best)
naive, prior, curves = naive_bayes(train, [0, 1, 2])
fixed, _, strength = naive_bayes(train, [3])
lr = logistic(train)
for name, model in (("base rates", base_rates(train)), ("naive Bayes, three clues", naive), ("naive Bayes, one strength clue", fixed), ("logistic regression", lr)):
    right = sum(max(RESULTS, key=model(r["x"]).get) == r["y"] for r in test)
    picks = Counter(max(RESULTS, key=model(r["x"]).get) for r in test)
    print(f"test: {name}: log loss {log_loss(model, test):.4f}, accuracy {right / len(test):.1%}, picks {dict(picks)}")
gap = [log(lr(r["x"])[r["y"]]) - log(fixed(r["x"])[r["y"]]) for r in test]
avg = sum(gap) / len(gap)
print(f"one strength clue minus logistic regression: {avg:+.4f}, luck margin +/- {2 * sqrt(sum((g - avg) ** 2 for g in gap) / (len(gap) - 1) / len(gap)):.4f}")

# 3. counting the same news twice: test matches where three clues make a home win 80% or more
sure = [r for r in test if naive(r["x"])["H"] >= 0.8]
print(f"three clues say home 80%+ in {len(sure)} test matches: they average {sum(naive(r['x'])['H'] for r in sure) / len(sure):.0%}, "
      f"one strength clue {sum(fixed(r['x'])['H'] for r in sure) / len(sure):.0%}, home wins happened {sum(r['y'] == 'H' for r in sure) / len(sure):.0%}")
even = [0, 0, 0, 0]
print("even game: three clues", {c: f"{p:.1%}" for c, p in naive(even).items()}, "| one clue", {c: f"{p:.1%}" for c, p in fixed(even).items()})

# 4. clue by clue, for the test match where the two disagree most about a home win
r = max(test, key=lambda r: abs(naive(r["x"])["H"] - fixed(r["x"])["H"]))
print(f"{r['match']} ({r['date']}), gaps " + ", ".join(f"{v:+.2f}" for v in r["x"][:3]) + f", strength {r['x'][3]:+.2f}")
for k in range(4):
    p = naive(r["x"], upto=k)
    print(f"  naive, after {k} clue{'s' if k != 1 else ''}: home {p['H']:.0%}, draw {p['D']:.0%}, away {p['A']:.0%}")
p = fixed(r["x"])
print(f"  one strength clue: home {p['H']:.0%}, draw {p['D']:.0%}, away {p['A']:.0%}")

# 5. what the model learned, for the Naive Bayes Match Predictor
print("base rates:", {c: round(prior[c] / len(train), 4) for c in RESULTS})
for j, label in enumerate(CLUES + ["strength"]):
    source = strength if j == 3 else curves
    print(f"{label}:", {c: tuple(round(v, 4) for v in source[c, j]) for c in RESULTS})