Premier League: team strengths, a title race and a fight with the bookmaker

#The brief

You have joined the analytics desk of a sports media company. The editor wants two things:

1. Team strength estimates - attack and defence ratings with honest uncertainty, not a league table with extra steps.
2. Probabilistic forecasts - "Arsenal have an X% chance of the title, West Ham a Y% chance of going down" - that the desk can publish mid-season without embarrassing itself.

The head of the desk is a sceptic: "Bookmakers do this for a living. Show me your numbers next to theirs on matches your model has not seen, then we talk."

You will build the model on the complete 2024/25 season, then deploy it on 2025/26 as it stood on 13 February 2026 - 26 rounds played, 12 to go - and forecast the rest. We know how the season ended, so we can mark your homework.

#How this notebook works

• Each task states what to deliver, not how. Write your code in the YOUR CODE HERE cells.
• Stuck? h.hint("task2") reveals hints one level at a time: nudge → approach → code skeleton. Try to get by on nudges.
• h.check("task2", home=...) compares your numbers with the reference solution.
• A full worked solution lives in notebooks/solutions/. Open it only when you are done (or truly stuck).

#The data

One row per match. The columns you need: HomeTeam, AwayTeam, full-time home and away goals FTHG and FTAG, the full-time result FTR (H, D or A), and Bet365's pre-match decimal odds B365H, B365D, B365A (a winning 1-unit bet at odds 2.5 pays back 2.5). There are many more columns (shots, cards, other bookmakers); ignore them or use them in "Going further".

A helper for the bookkeeping, so you can spend your time on the modelling. Ties on points are broken by goal difference, then goals scored.

#Task 0 · Look at the data

Use the 2024/25 season.

Deliver

1. How large is home advantage - in goals and in results?
2. Are goals plausibly Poisson? Compare the observed distribution of home goals and of away goals with a Poisson of the same mean, and compare means with variances.
3. One sentence: is a Poisson likelihood a defensible starting point?

Home teams scored 1.51 goals per match against 1.42 for away teams, and won 41% of matches against 35% - a real but, in this particular season, modest home advantage.

The marginal goal distributions sit close to a Poisson with the same mean, and variances are close to the means (1.63 vs 1.51 at home, almost identical away). Note that this is the marginal distribution over all fixtures: a mixture of Poissons with different rates (Liverpool at home to Southampton, and the reverse) must be somewhat overdispersed, so a variance slightly above the mean is what a Poisson model with team effects predicts. The same goes for the negative correlation between home and away goals (-0.14): strong teams score a lot and concede little, so it need not mean that goals within a match depend on each other. Poisson is a defensible starting point; we will check it properly in Task 3.

#Task 1 · The textbook model

The classic model for football scores gives every team an attack and a defence rating. For a match between home team $h$ and away team $a$:

$$ \begin{aligned} \text{home goals} &\sim \text{Poisson}(\lambda_h), & \log \lambda_h &= \mu + \eta + \text{att}_h - \text{def}_a \\ \text{away goals} &\sim \text{Poisson}(\lambda_a), & \log \lambda_a &= \mu + \text{att}_a - \text{def}_h \end{aligned} $$

with $\mu$ an intercept and $\eta$ the home advantage. A colleague suggests "uninformative priors, let the data speak": independent Normal(0, 3) on $\mu$, on $\eta$ and on every $\text{att}$ and $\text{def}$.

Deliver

1. Fit exactly this model to 2024/25 with default sampler settings, using coords/dims for teams.
2. Read the diagnostics. Something is wrong. Explain what, in terms of the model - not in terms of the sampler - and produce one plot or table that demonstrates your explanation.
3. Which quantities does this posterior pin down well in spite of the problem, and which not at all? Give one example of each, with numbers.

#Solution

No divergences - but r_hat of 1.05-1.08 and a bulk ESS of about 60 for the intercept and the team effects: four chains, 4000 draws, and the information of sixty. Yet home, in the same model, has an ESS of 2000 and r_hat of 1.00. That contrast is the clue.

Look at the posterior sd of the attack ratings: about 0.6 on the log scale, i.e. we supposedly cannot tell whether Liverpool score 0.5x or 2x what an average team does, after 38 matches. Now look at a difference between two teams:

The model is not identifiable. The likelihood only ever sees the sums $\mu + \text{att}_i - \text{def}_j$, and those are unchanged if you

• add a constant to every attack rating and subtract it from $\mu$, or
• add the same constant to every attack and every defence rating.

So there are two directions in parameter space along which the data say nothing at all. The posterior is a long, thin ridge (both panels): its length is set purely by the prior

• make the prior wider and the ridge gets longer - while its width is what the data determined. NUTS has to crawl along it, hence the dismal ESS. It would be a mistake to "fix" this with more draws or a higher target_accept: the sampler is reporting a real property of the model.

What is pinned down: differences between teams (sd about 0.22 instead of 0.6), the home advantage, and the overall scoring level $\mu + \overline{\text{att}} - \overline{\text{def}}$ (sd 0.04, against 0.85 for $\mu$ on its own). Those are the identified quantities, and the fix is to parameterise the model so that only they are free.

#Task 2 · An identified, hierarchical model

Deliver

1. A reformulated model in which every parameter is identified by construction (not by luck of the prior), and in which the spread of attacking and of defensive ability across the league is learned from the data rather than fixed by you.
2. Priors you can defend, backed by a prior predictive check on a scale a football fan would recognise.
3. Clean diagnostics, and a comparison with Task 1: what happened to the effective sample size?
4. A plot of attack and defence ratings with uncertainty. Then compare your attack ratings with the naive estimate "log of goals scored relative to the league average". Which way does the model move teams, and why?

#Solution

pm.ZeroSumNormal constrains the ratings to sum to zero across teams, which removes both flat directions at once: attack ratings are now "relative to the average team", and the intercept is unambiguously the log scoring rate of an average away team against an average defence. Giving the ratings a learned scale sigma_att, sigma_def makes the model hierarchical: teams are shrunk towards the league average by an amount the data choose.

Priors, on the log scale: teams score around 1-2 goals a match, so intercept ~ Normal(0, 0.5) covers rates from about 0.4 to 2.7. Home advantage multiplies scoring by maybe 1.1-1.4; Normal(0, 0.25). A rating sd of 0.5 would already mean the best attacks score $e^{1} \approx 2.7$ times the average, so HalfNormal(0.5) is generous.

The model is wrapped in a function with pm.Data inputs because we will fit it again to another season and predict fixtures it has not seen. (The noncentred switch is for Task 4.)

A median of one goal, nine in ten simulated performances at three goals or fewer, and just under 2% at more than seven. No side scored more than seven in either season of our data, so the prior is generous in the tail - which is what weakly informative means - but nothing absurd dominates. Fit:

Same likelihood, same data; ESS is 50 to 150 times higher for the intercept and the ratings, because the ridge is gone. (home was identified all along and barely changes.)

The home advantage in 2024/25 is small: about +6% goals, with an interval that includes zero. Team abilities have an sd of roughly 0.25 (attack) and 0.22 (defence) on the log scale.

Liverpool's attack and Arsenal's defence stand out at the top; the three relegated sides (Southampton, Ipswich, Leicester) are clearly the worst at both ends. But intervals are wide: after a full season, most mid-table teams cannot be separated.

The model's ratings are a compressed version of the raw ones, by about 30% (slope 0.69): extreme scoring records are partly luck, and the hierarchical prior says so. Southampton's 26 goals in 38 matches were bad, but the model's best guess is that their attack was not quite that bad. Nobody chose the 30%: it follows from sigma_att and from how much a season of goals can tell you about one team.

#Task 3 · Check the model against things it was not fitted to

The model was fitted to home goals and away goals, separately. Football is about results. A model can match both goal distributions and still get the things people care about wrong.

Deliver posterior predictive checks, each with a verdict, for

1. the share of draws in a season,
2. the frequency of the most common scorelines (0-0, 1-1, 1-0, ...),
3. one further statistic of your choice that the likelihood does not target directly.

Does anything fail badly enough to matter for forecasting results?

The rootograms cover what the model was fitted to: observed counts (black) sit inside or at the edge of the predictive intervals. Now for what it was not. Each posterior predictive draw is one replicated season of 380 matches:

Verdict. The share of draws - the quantity most likely to embarrass a model built from two independent Poissons - is reproduced: the model expects about 23%, the season had 24.5%, well inside the replicated range. Of the eleven most common scorelines, ten fall inside their 94% intervals; only 2-2 was clearly more frequent than the model expects, and 0-0 was rarer than expected, though not alarmingly so.

The one check that comes close to failing is the third: the variance of total goals per match was 2.6, and only about 2% of replicated seasons are that tame. The fourth panel says where it comes from. Team effects alone do make home and away goals negatively correlated in the replications (about -0.06), as we argued in Task 0, but the real season had -0.14, further out than roughly 94% of replications. When one side scored more, the other scored a little less than independence allows. Goals within a match are mildly dependent - the finding that made Dixon & Coles (1997) famous, although their version (too many 0-0 and 1-1) is not quite what this season shows.

Does it matter for forecasting results? The draw share and the scorelines say: not much. Two tail probabilities of 0.02 and 0.06 among several checks on one season are a lead, not a verdict. We proceed, and keep the dependence between the two scores in mind as the place to look if the forecasts disappoint.

#Task 4 · Forecast the run-in

Now the live season. It is 13 February 2026: every team has played 26 of 38 matches. played is what you may fit to; remaining holds the 120 fixtures still to come. (It also holds their results and odds. Those are for marking - do not let the model see them.)

Deliver

1. Your Task 2 model fitted to played, with clean diagnostics.
2. Simulations of the 120 remaining fixtures from the posterior predictive (parameter uncertainty and match-day randomness), turned into complete final league tables - one table per posterior draw, with the usual tie-breakers.
3. For every team: expected final points with a 94% interval, P(title), P(top 4), P(relegation) (bottom three).
4. The reckoning: compare with league_table(epl26). Were the favourites right? How many of the 20 final points totals fall inside your 94% intervals, and how many should?

#Solution

"Your Task 2 model, with clean diagnostics" hides a small trap. With 26 matches per team instead of 38, sigma_att and sigma_def are less well determined and values near zero become plausible - the funnel geometry of E02. Fit the model both ways and compare:

Neither fit diverges, but the centred one is visibly less comfortable: for the two sigmas a minimum bulk ESS of about 450 (tail: 320) and an r_hat of 1.008, against about 810 (610) and 1.004 for the non-centred version. With the full season of Task 2 it is the other way round (try it): the more the data say about every team, the better the centred form does. The centred fit would pass, but title and relegation probabilities lean on the tails of the posterior, so we take the better one.

One number stands out: home is 0.21, i.e. home sides score about 23% more - more than three times last season's estimate, with an interval that excludes zero. Whether home advantage really moves that much from one year to the next, or whether one of the two seasons is a fluke, is something a single-season model cannot know. We will see in Task 6 what the bookmaker thinks.

Now swap the fixtures in with pm.set_data and sample the posterior predictive for the unseen matches. We also keep the match rates mu_home, mu_away for Task 5.

Bookkeeping, vectorised over draws: start every simulated season from the real table on 13 February, add the simulated results, and rank with points, then goal difference, then goals scored (scaled so that one number sorts all three; a tiny random term breaks exact ties).

On 13 February the model made Arsenal clear but not overwhelming title favourites (about 70%, Manchester City most of the rest), and named Wolves, Burnley and West Ham as the most likely relegated sides, West Ham at roughly 60%. Arsenal won the league and exactly those three went down. The model's four most likely top-four finishers were the actual top four, but do not read much into that: the last place was close to a coin toss between Manchester United (44%) and Chelsea (41%), and it happened to fall the model's way.

More important than the hits: calibration. 18 of 20 final totals fall inside their 94% intervals, where we would expect about 19. The two misses are instructive. Chelsea took 8 points from their last 12 matches after 44 from 26, and Burnley took 4. At the other end, Manchester United's 26 points from 12 matches only just stayed inside their interval. The model assumes a team's strength is constant over the season; form, injuries, managerial changes and fixture congestion are all outside it.

Note how wide the intervals are: roughly ±9 points with 12 matches to go. Most of that is irreducible match-day randomness, not parameter uncertainty: fixing every team's strength at its posterior mean only reduces the sd of final points from 4.8 to 4.4. Better ratings would not make these intervals much narrower. A desk that publishes "Arsenal will finish on 79 points" rather than "between 70 and 88" is selling false precision.

#Task 5 · Your numbers next to the bookmaker's

The sceptic's test. For each of the 120 held-out matches, the model and the bookmaker each give three probabilities (home win, draw, away win), and then one of the three happened.

Two proper scoring rules, both "lower is better", averaged over matches, where $p_{ik}$ is the probability given to outcome $k$ in match $i$ and $o_{ik}$ is 1 if it happened, else 0:

$$\text{log score} = -\frac{1}{N}\sum_i \log p_{i,\text{observed}} \qquad\qquad \text{Brier} = \frac{1}{N}\sum_i \sum_{k \in \{H,D,A\}} (p_{ik} - o_{ik})^2$$

Deliver

1. Model probabilities for H/D/A for every held-out match, from the Task 4 posterior.
2. Bookmaker probabilities from the Bet365 odds. Raw inverse odds do not sum to one; how big is the bookmaker's margin (the "overround"), and what do you do about it?
3. A reference forecaster that knows nothing about teams.
4. A table of both scores for all three forecasters - and a statement about whether 120 matches are enough to separate them. Quantify it.

#Solution

The outcome probabilities could be read off the simulated scores of Task 4, but with 4000 draws they would carry Monte Carlo noise of up to about 0.008 each. We have the posterior draws of the scoring rates for each fixture, so we can do without that noise and compute the probabilities exactly: for each draw the score matrix is an outer product of two Poisson pmfs; average over draws (that is the posterior predictive), then sum the lower triangle (home win), diagonal (draw) and upper triangle (away win).

Bet365's odds imply probabilities that sum to about 1.055: a 5.5% margin. Dividing by the sum (proportional normalisation) is the simplest way to remove it.

The bookmaker wins, on both scores. Our model sits between the bookmaker and the know-nothing baseline, a little closer to the bookmaker. But look at the standard errors of the paired differences: each gap is about one to two standard errors. With 120 matches we cannot formally separate any pair of these forecasters - football results are so noisy that a useful edge in log score is a few hundredths, and measuring that takes many hundreds of matches. The honest summary for the head of desk: "as expected, no evidence that we beat the market; mild evidence that we beat nothing at all; ask me again in three seasons".

#Task 6 · Would you bet on it?

The head of desk reads your table and asks the natural follow-up: "Fine, the market is better on average. But when your model disagrees with it - is there money in that? And what does the market know that your model does not?"

Deliver

1. A back-test on the held-out matches: stake 1 unit on every outcome where your model sees positive expected value at the Bet365 odds (the real odds, margin included). Number of bets, profit, return on investment - with an uncertainty interval on the ROI.
2. A diagnosis of how the model's probabilities differ from the market's. Is the model more or less opinionated? If you mix the two forecasts, w * model + (1 - w) * market, which w scores best on the held-out matches, and what does that tell you?
3. A recommendation in three sentences. What is the model good for, what is it not good for, and what would you add first?

The back-test. The model "finds value" in 150 bets - which should already make you suspicious: a market this liquid does not misprice 40% of all outcomes. The realised ROI is a few percent negative, about what the bookmaker's margin predicts for someone betting without an edge, and the bootstrap interval runs from about -30% to +25%. For scale: blindly backing every outcome of every match lost only 1.4% on these matches, less than the margin - that is how noisy 120 matches are. They can neither prove nor disprove an edge; the prior expectation (no edge) stands.

The diagnosis. On average the two forecasters agree almost exactly: 44% home wins, 24% draws, 32% away wins. So the market, too, believes in this season's larger home advantage (and both expected fewer draws than the 29% that happened). The difference is in the spread, and the scatter is the most useful picture in this notebook. The two forecasters agree on the ordering of matches, but the cloud is flatter than the diagonal: where the market says 75%, the model says 60%; where the market says 15%, the model says 20-25%. The sd of the home-win probability across matches is 0.10 for the model and 0.17 for the market. Our model is less opinionated than the market, not more. That is the hierarchical prior at work: after 26 matches it shrinks Arsenal and Wolves towards the average, because from goals alone it cannot rule out luck. The market can, because it knows things our model does not: previous seasons, squad value, injuries, expected goals, line-ups. The best mixture puts a weight of only about 0.1 on the model - on these matches it adds next to nothing to the market's forecast.

Recommendation. Publish the model for what it is good at: transparent team ratings with honest uncertainty, and season simulations (title, top four, relegation) that the odds for single matches do not give you and that proved well calibrated. Do not use it to second-guess match odds, and do not bet with it. The diagnosis says our ratings are too timid, not too bold, so the first thing to add is information that sharpens them faster than goals can: shots or expected goals, and what we knew before the season started.

#Going further

• Last season as a prior. Fit both seasons jointly, with a 2025/26 rating centred on a fraction of the 2024/25 rating and a separate prior for promoted teams. Does it score better on the held-out matches? (A quick attempt of ours scored 1.056 against 1.052 for the single-season model at round 26 - no better - and its sampler diagnostics needed work. Does it help at round 8, when the current season says little?)
• Dependent scores. Task 3 found home and away goals more negatively correlated than the model allows, and the run-in had more draws than either forecaster expected. Add the low-score correction of Dixon & Coles (1997) with pm.Potential or pm.CustomDist, or a match-level random effect shared by both sides with opposite signs. Does the log score improve?
• Time-varying strength. Chelsea and Burnley broke the model in Task 4. Let ratings follow a pm.GaussianRandomWalk over match rounds. Is there evidence that form is real, or is it noise?
• Better data. Shots on target (HST, AST) are about three times as plentiful as goals. Model them jointly with goals and see whether ratings sharpen faster early in the season.