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

Source: https://www.footballdatascience.co.uk/learn/was-seven-nil-a-freak
Published: 2026-10-02

> Seven times since 2001/02 a Scottish Premiership side has won 7–0. A Bayesian model says each was a long shot, but also that a league should see more of them than it does. Checking the model against every score shows where it goes wrong, and how far a simple fix gets.

**On the terraces:** Seven-nil: a freak result, or just football? This piece asks a prediction model what scores it expects, finds it expects more thrashings than really happen, and shows why.

## The football question

On 31 March 2017 Aberdeen won **7–0 at Dundee**. Since 2001/02 a Scottish Premiership side has won 7–0 seven times: Celtic five times, Hibernian once and Aberdeen once. Bigger still were Celtic's two 9–0s, at home to Aberdeen in 2010 and away at Dundee United in 2022, and Rangers' 8–0 against Hamilton in 2020.

Were those **freaks**, scores a sensible model says shouldn't happen? Or just what a league throws up if you play enough matches?

## The concept

To answer, a model has to say more than who's likely to win. It has to say what **every** score should look like, across a whole league. That spread of possible scores is the model's **predictive distribution**, and comparing it with what really happened is a **posterior predictive check**: work out the scores your beliefs imply, then see whether reality looks like them.

The model is the one from [updating a team's scoring rate](/learn/updating-scoring-rates) and [simulating the league](/learn/simulating-the-league). Before every match, each side's attack and defence start from last season, counted as 20 matches, and are updated with every match so far this season. A side's expected goals are its attack times the other side's defence, adjusted for home or away, and the goals themselves are [Poisson](/learn/poisson-distribution). Run before every match from 2001/02 to 2025/26, that's **11,302 side-matches**: one for each team in each game.

There are two ways to turn the beliefs into scores:

- **The model:** use the best guess of each side's strength.
- **The model allowing for doubt:** we're never sure of those strengths, so draw plausible ones from the range the evidence allows, and average. This is the textbook Bayesian forecast. It spreads the scores out, with more blanks and more big wins; it's where the [Negative Binomial](/learn/negative-binomial-distribution) comes from.

I expected allowing for doubt to explain the thrashings. It doesn't.

## A football example

Every side-match from 2001/02 to 2025/26, how often a side scored five or more, and seven or more:

| 2001/02 to 2025/26 | 5 or more | 7 or more |
|---|---|---|
| Real | 215 | 18 |
| The model | 278 | 33 |
| Allowing for doubt | 309 | 43 |

The model expects **almost twice as many sevens as really happened**, and allowing for doubt makes it worse, 43. Blanks go the same way: 3,240 real, 3,396 from the model, 3,498 allowing for doubt. **The model is too spread out**, and the textbook fix spreads it further.

### What's wrong

It isn't the average. Over the 25 seasons sides scored 15,145 goals and the model expected 15,166. So the problem is the spread. Group every side-match by how many goals the model expected:

<figure class="rank-chart">
<div role="img" aria-label="Goals a game, model against real, by the model's expected goals. 0 to 0.5: model 0.42, real 0.54. 0.5 to 1: 0.80, 0.89. 1 to 1.5: 1.24, 1.25. 1.5 to 2: 1.70, 1.62. 2 to 2.5: 2.23, 2.22. 2.5 to 3: 2.72, 2.48. 3 or more: 3.36, 2.78.">

</div>
<figcaption>In the middle the model is spot on. At the edges it's too extreme: sides it expects to struggle score more than it thinks, and sides it expects to run riot score less.</figcaption>
</figure>

Where the model expects a side to score **3.36 a game, they score 2.78**. It's too sure about mismatches: multiply a strong attack by a weak defence and the answer overshoots. Thrashings come from mismatches, so it expects too many of them.

Do sides that are well ahead simply ease off? Not on this evidence. The side ahead at half-time scored more in the second half the bigger its lead: 0.72 goals level, 0.82 one up, 0.97 two up, 1.03 three up and 1.34 four or more up (only 38 matches). Big leads mostly belong to strong sides, so this doesn't prove nobody ever eases off, but it isn't the story here.

### A fix, and a fair test

If the model is too extreme, pull every forecast part of the way back towards the average:

$$\lambda' = c + b\,(\lambda - c)$$

<div class="plain" markdown="1">
In plain football

- **λ** is the model's expected goals for one side in one match.
- **c** is the average for a home side, or an away side, last season.
- **b** is how much of the gap from that average to keep: 1 trusts the model completely, 0 says every side is average.
- **λ′** is the adjusted expected goals.
</div>

It's the same idea as [partial pooling](/learn/partial-pooling): trust the extremes a little less. Following the rule for [choosing settings honestly](/learn/elo-tuned), **b** was chosen on the held-back seasons, 2016/17 to 2020/21, and came out at **0.88**: keep 88% of each forecast's distance from the average. Then the test seasons, 2021/22 to 2025/26, were checked once:

| Test seasons | Before | After |
|---|---|---|
| No goals (real 629) | 687 | 668 |
| 5 or more (real 47) | 62.1 | 55.4 |
| 7 or more (real 3) | 8.2 | 6.4 |

Every count moves towards reality, and the fit to the test seasons improves (log-likelihood −1.4460 before, −1.4413 after; higher is better). But it **closes less than half the gap**: the shrunk model still expects 55 scores of five or more where there were 47. A simple pull towards the average helps; it doesn't cure the overshoot.

## So, was 7–0 a freak?

The shrunk model's chance, before kick-off, that each 7–0 winner would score seven or more:

| 7–0 | Expected | Chance |
|---|---|---|
| Celtic 7–0 Aberdeen, 2002 | 2.67 | 1&nbsp;in&nbsp;51 |
| Hibernian 7–0 Livingston, 2006 | 2.52 | 1&nbsp;in&nbsp;67 |
| Celtic 7–0 St Mirren, 2009 | 2.28 | 1&nbsp;in&nbsp;111 |
| Celtic 7–0 Motherwell, 2016 | 2.92 | 1&nbsp;in&nbsp;34 |
| Dundee 0–7 Aberdeen, 2017 | 1.64 | 1&nbsp;in&nbsp;656 |
| Celtic 7–0 St Johnstone, 2019 | 2.24 | 1&nbsp;in&nbsp;123 |
| Celtic 7–0 St Johnstone, 2022 | 2.08 | 1&nbsp;in&nbsp;178 |

Expected: the winning side's expected goals in that match.

Every one was a long shot. But a league plays a lot of matches: over 25 seasons the shrunk model expects **25.9** scores of seven or more, **about one a season**. There were **18**. That's fewer, but luck alone would give 18 or fewer 6.6% of the time, so it isn't clear evidence that the model is still wrong about sevens.

**The verdict:** no 7–0 is a freak on its own, Aberdeen's at Dundee coming closest at 1 in 656. As a group they're what a league should expect, and if anything a little rarer.

## Why it matters

- **Check the whole spread, not just the winner.** A model can pick results well and still get the shape of the scores wrong, and score markets like over and under 2.5 goals live on that shape.
- **The textbook step isn't automatically right.** Allowing for doubt made things worse, because the trouble wasn't uncertainty about each side's strength but how two strengths combine.
- **Let the check point to the fix,** choose its setting on held-back seasons, and test it once.
- **"Freak" needs a number.** One match at 1 in 656 is rare; one a season across a league is not.

## Limitations

- **Goals only.** The model knows nothing of red cards, injuries or rotation, which make some thrashings more likely.
- **The fix is simple.** One shrink factor for every match; a model that combined attack and defence differently, or [Dixon-Coles](/learn/dixon-coles-ratings), might fit the edges better. Neither was tried here.
- **Two of the seven 7–0s fall in the held-back seasons** used to choose 0.88, so their chances in the table aren't fully out of sample.
- **The easing-off check is rough:** sides four up at half-time are few (38) and mostly strong.
- **Sevens are rare,** 18 in 25 seasons, so their own count is a small sample: the 6.6% above is that luck at work.

## Try it yourself

Take your team's biggest win this season and its usual goals a game at home. Put that into the [Poisson Match Predictor](/models/poisson) and read the chance of that many goals or more from the grid. Freak, or about once a season somewhere in the league?

## Reproduce the analysis

Download the Scottish Premiership files (SC0) for 2000/01 to 2025/26 from [football-data.co.uk](https://www.football-data.co.uk/scotlandm.php), saved as `SC0_0001.csv` and so on; they aren't rehosted on this site. It uses a fixed random seed and takes a few seconds. Then:

```python
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%}")
```

## Further reading

- [Posterior predictive distribution](https://en.wikipedia.org/wiki/Posterior_predictive_distribution), Wikipedia. Forecasting new data while carrying the uncertainty in the parameters through, the "allowing for doubt" version here.
- [Overdispersion](https://en.wikipedia.org/wiki/Overdispersion), Wikipedia. When data are more spread out than a model like Poisson allows; here the model was more spread out than the data.
- [Bayesian Data Analysis](https://sites.stat.columbia.edu/gelman/book/), by Gelman, Carlin, Stern, Dunson, Vehtari and Rubin. The authors' page for the book, with the whole book as a PDF for non-commercial use; its chapter on model checking sets out posterior predictive checks.