C02 · Premier League forecasting — Bayesian Modelling in Practice

The brief

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

Team strength estimates - attack and defence ratings with honest uncertainty, not a league table with extra steps.
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).

import arviz as az
import matplotlib.pyplot as plt
import numpy as np
import pandas as pd
import pymc as pm
from scipy import stats

from pymc_challenges import Hints, data

RANDOM_SEED = 2526
rng = np.random.default_rng(RANDOM_SEED)
az.style.use("arviz-variat")

h = Hints("C02")
h.tasks()

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

task0 - Look at the data (2 hints, 2 check(s))
task1 - The textbook model and what is wrong with it (3 hints, 1 check(s))
task2 - An identified, hierarchical model (3 hints, 3 check(s))
task3 - Posterior predictive checks on results (3 hints, 1 check(s))
task4 - Forecast the run-in (3 hints, 3 check(s))
task5 - Scoring against the bookmaker (3 hints, 3 check(s))
task6 - Would you bet on it? (3 hints, 1 check(s))

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".

COLUMNS = ["Date", "HomeTeam", "AwayTeam", "FTHG", "FTAG", "FTR", "B365H", "B365D", "B365A"]


def prepare(df):
    df = df[COLUMNS].copy()
    df["Date"] = pd.to_datetime(df.Date, dayfirst=True)
    return df.sort_values("Date", kind="stable").reset_index(drop=True)


data.describe("epl_2425")
epl25 = prepare(data.load("epl_2425"))
epl25.head()

epl_2425: All 380 match results with goals, shots and bookmaker odds.
  source: football-data.co.uk, English Premier League 2024/25.
  url:    https://www.football-data.co.uk/mmz4281/2425/E0.csv

	Date	HomeTeam	AwayTeam	FTHG	FTAG	FTR	B365H	B365D	B365A
0	2024-08-16	Man United	Fulham	1	0	H	1.60	4.20	5.25
1	2024-08-17	Ipswich	Liverpool	0	2	A	8.50	5.50	1.33
2	2024-08-17	Arsenal	Wolves	2	0	H	1.18	7.50	13.00
3	2024-08-17	Everton	Brighton	0	3	A	2.63	3.30	2.63
4	2024-08-17	Newcastle	Southampton	1	0	H	1.36	5.25	8.00

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.

def league_table(df):
    """League table (points, goals for/against, goal difference) from a frame of matches."""
    home = pd.DataFrame({"team": df.HomeTeam, "gf": df.FTHG, "ga": df.FTAG})
    away = pd.DataFrame({"team": df.AwayTeam, "gf": df.FTAG, "ga": df.FTHG})
    long = pd.concat([home, away])
    long["played"] = 1
    long["pts"] = np.select([long.gf > long.ga, long.gf == long.ga], [3, 1], 0)
    table = long.groupby("team")[["played", "pts", "gf", "ga"]].sum()
    table["gd"] = table.gf - table.ga
    table = table.sort_values(["pts", "gd", "gf"], ascending=False)
    table.insert(0, "pos", range(1, len(table) + 1))
    return table


league_table(epl25)

	pos	played	pts	gf	ga	gd
team
Liverpool	1	38	84	86	41	45
Arsenal	2	38	74	69	34	35
Man City	3	38	71	72	44	28
Chelsea	4	38	69	64	43	21
Newcastle	5	38	66	68	47	21
Aston Villa	6	38	66	58	51	7
Nott'm Forest	7	38	65	58	46	12
Brighton	8	38	61	66	59	7
Bournemouth	9	38	56	58	46	12
Brentford	10	38	56	66	57	9
Fulham	11	38	54	54	54	0
Crystal Palace	12	38	53	51	51	0
Everton	13	38	48	42	44	-2
West Ham	14	38	43	46	62	-16
Man United	15	38	42	44	54	-10
Wolves	16	38	42	54	69	-15
Tottenham	17	38	38	64	65	-1
Leicester	18	38	25	33	80	-47
Ipswich	19	38	22	36	82	-46
Southampton	20	38	12	26	86	-60

Task 0 · Look at the data

Use the 2024/25 season.

Deliver

How large is home advantage - in goals and in results?
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.
One sentence: is a Poisson likelihood a defensible starting point?

Your turn Work this task in the notebook.

Hints for this task — open one at a time

1 Nudge

Home advantage shows up twice: in goals (FTHG against FTAG) and in results (FTR). For the Poisson question, remember the one property that defines a Poisson: what is the relationship between its mean and its variance?

2 Approach

epl25.FTR.value_counts(normalize=True) gives the result shares. scipy.stats.poisson.pmf(np.arange(8), mean) gives the Poisson probabilities to draw on top of the observed shares epl25.FTHG.value_counts(normalize=True).

Think before you judge the variance: all 380 matches do not share one scoring rate. A mixture of Poissons with different rates has a variance above its mean. So mild overdispersion in the pooled goals is not evidence against a Poisson model with team effects.

Self-checks in the notebook (h.check("task0", …)):

home_goals_mean — mean home goals per match in 2024/25, about 1.51
away_goals_mean — mean away goals per match in 2024/25, about 1.42

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

Fit exactly this model to 2024/25 with default sampler settings, using coords/dims for teams.
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.
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.

Your turn Work this task in the notebook.

Hints for this task — open one at a time

1 Nudge

Compare the rows of az.summary. Some parameters mix badly, at least one mixes well. What distinguishes them? Then a thought experiment: add 5 to every attack rating and subtract 5 from the intercept. What happens to every lambda in the model - and so to the likelihood?

2 Approach

The likelihood only ever sees intercept + attack[i] - defence[j]. Shifting all attacks by c and the intercept by -c leaves it unchanged; so does shifting all attacks and all defences by c. Along those two directions only the Normal(0, 3) prior holds the posterior in place: a long flat ridge that NUTS has to crawl along.

Demonstrate it with the posterior draws (naive_idata.posterior is an xarray Dataset): a scatter of one team's attack rating against another's, or of the intercept against attack.mean("team") - defence.mean("team"). Then compare the posterior sd of a single rating with the sd of a difference between two ratings. Differences, the home advantage and intercept + mean(attack) - mean(defence) are what the data identify.

Teams: pd.Categorical(df.HomeTeam, categories=teams).codes gives integer indices that are consistent between the home and away columns.

3 Skeleton

teams = np.sort(epl25.HomeTeam.unique())
home_idx = pd.Categorical(epl25.HomeTeam, categories=teams).codes
away_idx = pd.Categorical(epl25.AwayTeam, categories=teams).codes

with pm.Model(coords={"team": teams}) as naive_model:
    intercept = pm.Normal("intercept", 0, 3)
    home = pm.Normal("home", 0, 3)
    attack = pm.Normal("attack", 0, 3, dims="team")
    defence = pm.Normal("defence", 0, 3, dims="team")
    pm.Poisson("home_goals", pm.math.exp(intercept + home + attack[___] - defence[___]), observed=___)
    pm.Poisson("away_goals", pm.math.exp(___), observed=___)
    naive_idata = pm.sample(random_seed=RANDOM_SEED)

post = naive_idata.posterior
diff = post["attack"].sel(team="Liverpool") - post["attack"].sel(team="Southampton")
print(post["attack"].sel(team="Liverpool").std(), diff.std())

Self-checks in the notebook (h.check("task1", …)):

sd_attack_difference — posterior sd of attack[Liverpool] - attack[Southampton] in the NAIVE model, about 0.22 (far smaller than the sd of either rating on its own)

Task 2 · An identified, hierarchical model

Deliver

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.
Priors you can defend, backed by a prior predictive check on a scale a football fan would recognise.
Clean diagnostics, and a comparison with Task 1: what happened to the effective sample size?
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?

Your turn Work this task in the notebook.

Hints for this task — open one at a time

1 Nudge

Two flat directions need two constraints. What is the most natural reference point for "attack strength" - a particular team, or the average team? And who should decide how different the teams of a league are: you, with a fixed prior sd, or the 380 matches?

2 Approach

pm.ZeroSumNormal("attack", sigma=sigma_att, dims="team") is a Normal constrained to sum to zero over team. One for attack and one for defence removes both flat directions; the intercept becomes the log scoring rate of an average away side. Putting a prior on sigma_att and sigma_def (e.g. HalfNormal) makes the model hierarchical.

Prior scales live on the log-rate: exp(0.5) is already 1.65x. For the prior predictive check look at simulated goals per team per match (idata.prior_predictive["home_goals"]): how often does a side score 8 or more?

You will refit this model on another season and predict unseen fixtures, so wrap it in a function, register the team indices with pm.Data(..., dims="match") and give the likelihoods shape=mu.shape. For the comparison with Task 1, the ess_bulk column of az.summary for both fits says it all.

3 Skeleton

def build_model(df, teams):
    coords = {"team": teams, "match": df.index.values}
    with pm.Model(coords=coords) as model:
        home_idx = pm.Data("home_idx", pd.Categorical(df.HomeTeam, categories=teams).codes, dims="match")
        away_idx = pm.Data("away_idx", ___, dims="match")
        intercept = pm.Normal("intercept", 0, ___)
        home = pm.Normal("home", 0, ___)
        sigma_att = pm.HalfNormal("sigma_att", ___)
        sigma_def = pm.HalfNormal("sigma_def", ___)
        attack = pm.ZeroSumNormal("attack", sigma=___, dims="team")
        defence = pm.ZeroSumNormal("defence", sigma=___, dims="team")
        mu_home = pm.Deterministic("mu_home", pm.math.exp(___), dims="match")
        mu_away = pm.Deterministic("mu_away", pm.math.exp(___), dims="match")
        pm.Poisson("home_goals", mu_home, observed=df.FTHG.values, dims="match", shape=mu_home.shape)
        pm.Poisson("away_goals", mu_away, observed=df.FTAG.values, dims="match", shape=mu_away.shape)
    return model

Raw attack rating to compare with: np.log(goals_for / goals_for.mean()) from league_table.

Self-checks in the notebook (h.check("task2", …)):

home — posterior mean of the home advantage on the LOG scale for 2024/25, about 0.06 (small that season)
sigma_att — posterior mean of the sd of attack ratings across teams (log scale), about 0.25
sigma_def — posterior mean of the sd of defence ratings across teams (log scale), about 0.22

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

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

Does anything fail badly enough to matter for forecasting results?

Your turn Work this task in the notebook.

Hints for this task — open one at a time

1 Nudge

One posterior predictive draw is a complete replicated season: 380 home scores and 380 away scores. Any statistic you can compute for the real season you can compute for each replicated one. Where does the real value sit among the replicated ones?

2 Approach

pm.sample_posterior_predictive(idata, extend_inferencedata=True), then az.extract(idata, group="posterior_predictive", var_names=["home_goals", "away_goals"]) gives arrays of shape (match, sample). The share of draws in replicated season s is (hg[:, s] == ag[:, s]).mean() - vectorise with .mean(axis=0).

az.plot_ppc_rootogram checks the two goal distributions the model was fitted to; the draws and scorelines need the joint behaviour of home and away goals, so compute those by hand. Summarise each check as observed value, replicated mean and a tail probability P(rep >= obs). Candidates for the third statistic: share of 0-0, variance of total goals, correlation of home and away goals, the biggest winning margin of the season.

3 Skeleton

pp = az.extract(idata25, group="posterior_predictive", var_names=["home_goals", "away_goals"])
hg, ag = pp["home_goals"].values, pp["away_goals"].values          # (match, sample)

draws_rep = (___).mean(axis=0)                                     # one number per replicated season
draws_obs = (epl25.FTHG == epl25.FTAG).mean()
print(draws_obs, draws_rep.mean(), (draws_rep >= draws_obs).mean())

share_1_1 = ((hg == 1) & (ag == 1)).mean(axis=0)                   # repeat for each scoreline

Self-checks in the notebook (h.check("task3", …)):

draw_share_replicated — mean share of drawn matches over REPLICATED 2024/25 seasons, about 0.23 (the observed share was 0.245)

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.)

data.describe("epl_2526")
epl26 = prepare(data.load("epl_2526"))

CUTOFF = pd.Timestamp("2026-02-13")
played = epl26[epl26.Date < CUTOFF].reset_index(drop=True)
remaining = epl26[epl26.Date >= CUTOFF].reset_index(drop=True)

print(f"{len(played)} matches played, {len(remaining)} remaining")
league_table(played)

epl_2526: All 380 match results with goals, shots and bookmaker odds.
  source: football-data.co.uk, English Premier League 2025/26.
  url:    https://www.football-data.co.uk/mmz4281/2526/E0.csv
260 matches played, 120 remaining

	pos	played	pts	gf	ga	gd
team
Arsenal	1	26	57	50	18	32
Man City	2	26	53	54	24	30
Aston Villa	3	26	50	37	27	10
Man United	4	26	45	47	37	10
Chelsea	5	26	44	47	30	17
Liverpool	6	26	42	41	35	6
Brentford	7	26	40	40	35	5
Everton	8	26	37	29	30	-1
Bournemouth	9	26	37	43	45	-2
Newcastle	10	26	36	37	37	0
Sunderland	11	26	36	27	30	-3
Fulham	12	26	34	35	40	-5
Crystal Palace	13	26	32	28	32	-4
Brighton	14	26	31	34	34	0
Leeds	15	26	30	36	45	-9
Tottenham	16	26	29	36	37	-1
Nott'm Forest	17	26	27	25	38	-13
West Ham	18	26	24	32	49	-17
Burnley	19	26	18	28	51	-23
Wolves	20	26	9	16	48	-32

Deliver

Your Task 2 model fitted to played, with clean diagnostics.
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.
For every team: expected final points with a 94% interval, P(title), P(top 4), P(relegation) (bottom three).
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?

Your turn Work this task in the notebook.

Hints for this task — open one at a time

1 Nudge

A forecast of the final table needs two kinds of uncertainty: you do not know the team strengths exactly, and even if you did, football is random. Which PyMC function gives you both at once for matches the model has not seen? Each draw it returns is one possible future - finish the season in each.

2 Approach

Build the model on played (all 20 teams of 2025/26 as coords), sample, then pm.set_data({"home_idx": ..., "away_idx": ...}, coords={"match": remaining.index.values}) and pm.sample_posterior_predictive(idata, predictions=True, var_names=["home_goals", "away_goals", "mu_home", "mu_away"]). Keep the rates - Task 5 can use them. az.extract(forecast, group="predictions") gives (match, sample) arrays.

Look at ESS and r_hat of sigma_att and sigma_def before you move on. Fewer matches mean the sigmas are less well determined, and that is when the funnel of E02 shows up. The non-centred form of a zero-sum effect is sigma * pm.ZeroSumNormal("z", sigma=1, dims="team") wrapped in a pm.Deterministic.

For the tables: start from league_table(played) (points, gd, gf per team), loop over the 120 fixtures and add points and goals for all draws at once with NumPy. Rank within each draw by points, then goal difference, then goals scored; a double argsort along the team axis turns a sort key into ranks. P(title) is (rank == 1).mean(axis=1).

3 Skeleton

pts = np.repeat(start.pts.values[:, None], n_draws, axis=1).astype(float)   # (team, sample)
for k in range(len(remaining)):
    h_k, a_k = home_codes[k], away_codes[k]
    pts[h_k] += np.select([hg[k] > ag[k], hg[k] == ag[k]], [3, 1], 0)
    pts[a_k] += ___
    # same idea for goal difference and goals scored
key = pts * 1e6 + (gd + 500) * 1e3 + gf
rank = (-key).argsort(axis=0).argsort(axis=0) + 1                          # 1 = champions
p_title = (rank == 1).mean(axis=1)
p_relegation = (___).mean(axis=1)

Self-checks in the notebook (h.check("task4", …)):

p_title_arsenal — probability on 13 February 2026 that Arsenal win the league (a probability between 0 and 1; about 0.7 - clear favourites, far from certain)
p_relegation_west_ham — probability on 13 February 2026 that West Ham finish in the bottom three, about 0.6
expected_points_arsenal — Arsenal's expected FINAL points total (57 already won plus the simulated run-in), about 79

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

Model probabilities for H/D/A for every held-out match, from the Task 4 posterior.
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?
A reference forecaster that knows nothing about teams.
A table of both scores for all three forecasters - and a statement about whether 120 matches are enough to separate them. Quantify it.

Your turn Work this task in the notebook.

Hints for this task — open one at a time

1 Nudge

Decimal odds of 2.5 mean the bookmaker pays as if the probability were 1/2.5. Add up the three implied probabilities of a match - why is the sum above one, and who keeps the difference? And for the "how sure are you" part: both forecasters are scored on the same 120 matches, so what is the right way to compare them?

2 Approach

Bookmaker: q = 1 / odds, overround q.sum(axis=1) - 1, then divide each row by its sum (proportional normalisation; fancier methods exist, the differences are small).

Model: P(H), P(D), P(A) can be counted from the simulated scores of Task 4, or computed with less Monte Carlo noise from the rate draws: for each draw the score matrix is the outer product of two Poisson pmfs; average over draws, then sum below / on / above the diagonal.

Reference: the H/D/A frequencies in played, the same three numbers for every match.

Uncertainty: compute per-match score differences between two forecasters; their mean divided by std / sqrt(N) tells you how many standard errors apart they are.

3 Skeleton

g = np.arange(13)
pmf_h = stats.poisson.pmf(g[None, :, None], mu_home[:, None, :])     # (match, goals, sample)
pmf_a = stats.poisson.pmf(g[None, :, None], mu_away[:, None, :])
score = np.einsum("mis,mjs->mij", pmf_h, pmf_a) / mu_home.shape[1]   # (match, home goals, away goals)
p_home = np.tril(score, -1).sum((1, 2))
p_draw = ___
p_away = ___

q = 1 / remaining[["B365H", "B365D", "B365A"]].values
p_book = q / ___

outcome = remaining.FTR.map({"H": 0, "D": 1, "A": 2}).values
log_score = -np.log(p[np.arange(len(outcome)), outcome]).mean()

Self-checks in the notebook (h.check("task5", …)):

overround — mean Bet365 overround on the held-out matches as a fraction (sum of inverse odds minus 1), about 0.055
log_score_bookmaker — mean log score (negative log probability of the observed result, natural log) of the normalised Bet365 probabilities on the 120 held-out matches, about 1.035
log_score_model — mean log score of your model on the same matches, about 1.05 - a little above (worse than) the bookmaker's and below log(3) = 1.0986

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

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.
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?
A recommendation in three sentences. What is the model good for, what is it not good for, and what would you add first?

Your turn Work this task in the notebook.

Hints for this task — open one at a time

1 Nudge

A bet at decimal odds o on an event you believe has probability p has expected profit p * o - 1 per unit staked. How many such "value bets" does your model find? Is that number plausible for one of the most liquid betting markets in the world? Then plot model probability against market probability and ask which of the two is bolder.

2 Approach

Profit per bet: odds - 1 if the outcome happened, -1 if not. ROI is total profit over total stakes. For the interval, bootstrap over matches (resample the 120 rows with replacement), not over individual bets - bets on the same match are dependent.

For the diagnosis: scatter p_model against p_bookmaker. A cloud flatter than the diagonal means the model is less opinionated than the market (shrinkage!). For the mixture, evaluate the held-out log score of w * p_model + (1 - w) * p_bookmaker on a grid of w. A best w near zero says the model adds nothing the market does not know already. (This w is tuned on the test matches, so it is a diagnostic, not a strategy.)

3 Skeleton

odds = remaining[["B365H", "B365D", "B365A"]].values
value = p_model * odds - 1 > 0                                   # (match, outcome) bets placed
profit = np.where(onehot == 1, ___, ___) * value

boot = []
for _ in range(5000):
    i = rng.integers(0, len(remaining), len(remaining))
    boot.append(profit[i].sum() / value[i].sum())

weights = np.linspace(0, 1, 41)
blend = [log_score(w * p_model + (1 - w) * p_bookmaker) for w in weights]
best_w = weights[np.argmin(blend)]

Self-checks in the notebook (h.check("task6", …)):

best_blend_weight — weight on the model in the model/market mixture with the best held-out log score - small (about 0.1), the market dominates

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.

h.progress()

C02: 0/20 hints used, 14/14 checks passed.

Reference solution

A complete, executed solution with commentary on every task — including the ways the obvious model fails. Try the hint ladder first.

Open the worked solution →

#The brief

#How this notebook works

#The data

#Task 0 · Look at the data

#Task 1 · The textbook model

#Task 2 · An identified, hierarchical model

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

#Task 4 · Forecast the run-in

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

#Task 6 · Would you bet on it?

#Going further

Reference solution

The brief

How this notebook works

The data

Task 0 · Look at the data

Task 1 · The textbook model

Task 2 · An identified, hierarchical model

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

Task 4 · Forecast the run-in

Task 5 · Your numbers next to the bookmaker's

Task 6 · Would you bet on it?

Going further